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Abstract 


Naturally  occurring  radio  noise  above  approximately  100  MHz  is  well  modeled  for 
most  applications  as  a  Gaussian  random  process;  however,  atmospheric  radio  noise 
below  100  MHz  is  often  impulsive  in  nature  and  is  not  well  modeled  as  Gaussian. 
Atmospheric  events  (mainly  lightning  strokes,  which  create  electromagnetic  emis¬ 
sions  known  as  sferics)  produce  large,  clustered  impulses  in  the  noise  waveform  seen 
at  a  receiving  antenna,  causing  the  waveform  to  vary  greatly  from  typical  Gaus¬ 
sian  background  noise.  Due  to  large  variations  in  sferic  activity  on  a  seasonal  and 
diurnal  basis  and  with  the  passing  of  individual  storms,  atmospheric  noise  is  also 
non-stat  ionary. 

The  objective  of  this  work  is  the  statistical  characterization  and  modeling  of 
atmospheric  radio  noise  in  the  range  10  Hz  -  60  kHz  (denoted  low-frequency  noise), 
with  the  specific  goal  of  improving  communication  systems  operating  in  this  range. 
The  analyses  are  based  on  many  thousands  of  hours  of  measurements  made  by  the 
Stanford  Radio  Noise  Survey  System. 

The  statistics  analyzed  include  seasonal  and  diurnal  variations,  amplitude  prob¬ 
ability  distributions  (APDs),  impulse  interarrival  time  distributions,  background 
noise  statistics,  and  power  spectral  densities.  Noise  characteristics  derived  from 
these  analyses  are  presented,  and  several  noise  models  that  accurately  represent 
low-frequency  noise  APDs  are  compared.  The  most  accurate  model  for  representing 
low-frequency  noise  APDs  is  found  to  depend  on  geographic  location,  time  of  year 
and  day,  bandwidth,  and  center  frequency,  but  two  of  the  simplest  models  (i.e., 
each  with  only  two  parameters)  are  found  to  give  extremely  good  performance  in 
general.  These  are  the  Hall  and  alpha-stable  (or  a-stable)  models,  both  of  which 
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approximate  the  Rayleigh  distribution  for  low  amplitude  values  but  decay  with  an 
inverse  power  law  for  high  amplitude  values.  It  is  concluded  that  the  Hall  model 
is  the  optimal  choice  in  terms  of  accuracy  and  simplicity  for  locations  exposed  to 
heavy  sferic  activity  ( e.g .,  low-latitude  regions),  and  the  a-stable  model  is  best  for 
locations  relatively  distant  from  heavy  sferic  activity  (e.g.,  high-latitude  regions). 

Based  on  the  statistical  characteristics  of  the  noise  data,  a  new  clustering  Poisson 
atmospheric  noise  model  is  developed.  This  model  is  based  on  several  previously 
known  statistical-physical  models,  but  in  addition  takes  into  account  the  clustering 
of  sferic  impulses.  It  is  shown  that  the  clustering  model  accurately  predicts  the 
statistical  features  found  in  low-frequency  radio  noise  data. 

Finally,  since  many  low-frequency  digital  communication  systems  operate  in 
high-latitude  regions,  which  are  relatively  distant  from  heavy  sferic  activity,  it  is 
of  value  to  compare  the  bit  error  rate  (BER)  performance  of  receivers  specifically 
designed  for  a-stable  noise  with  the  BER  performance  of  conventional  low-frequency 
receivers.  The  communication  signal  formats  examined  are  quadrature  phase-shift 
keying  (QPSK)  and  16  point  quadrature  amplitude  modulation  (16  QAM).  Hun¬ 
dreds  of  simulations  using  time-series  data  from  various  times  and  locations  and  at 
various  center  frequencies  and  bandwidths  are  performed,  and  the  following  results 
are  found  uniformly:  for  QPSK  signals,  virtually  no  performance  improvement  is 
gained  when  using  an  a-stable  receiver  instead  of  the  best  conventional  receiver,  but 
for  16  QAM  signals,  an  improvement  of  several  dB  is  gained  by  using  an  a-stable 
receiver. 
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Chapter  1 
Introduction 


Naturally  occurring  radio  noise  above  approximately  100  MHz  is  well  modeled  for 
most  applications  as  a  Gaussian  random  process;  however,  atmospheric  radio  noise 
below  100  MHz  is  often  impulsive  in  nature  and  is  not  well  modeled  as  Gaussian. 
Atmospheric  events  (mainly  lightning  strokes,  which  create  electromagnetic  emis¬ 
sions  known  as  sferics )  produce  large,  clustered  impulses  in  the  noise  waveform  seen 
at  a  receiving  antenna,  causing  the  waveform  to  vary  greatly  from  typical  Gaus¬ 
sian  background  noise.  Due  to  large  variations  in  sferic  activity  on  a  seasonal  and 
diurnal  basis  and  with  the  passing  of  individual  storms,  atmospheric  noise  is  also 
non-stationary. 

The  objective  of  this  work  is  the  statistical  analysis,  characterization  and  mod¬ 
eling  of  atmospheric  radio  noise  in  the  frequency  range  10  Hz  -  60  kHz,  with  the 
specific  goal  of  improving  communication  systems  operating  in  this  range.  The  anal¬ 
yses  are  based  on  many  thousands  of  hours  of  measurements  made  by  the  Stanford 
Radio  Noise  Survey  System. 

This  chapter  begins  with  an  overview  of  10  Hz  -  60  kHz  radio  noise,  which 
includes  most  or  all  of  the  extremely-low  frequency  (ELF;  3  Hz  -  3  kHz),  very-low 
frequency  (VLF;  3  kHz  -  30  kHz)  and  low  frequency  (LF;  30  kHz  -  300  kHz)  bands. 
Following  this  overview  is  a  description  of  the  measurement  systems  used  to  collect 
the  noise  data.  The  contents  of  each  chapter  axe  then  summarized,  and  finally, 
the  main  contributions  of  the  dissertation  are  presented.  From  this  point  on,  the 
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terms  ELF,  VLF  and  LF  refer  to  the  individual  frequency  bands,  while  the  term 
low-frequency  refers  specifically  to  the  10  Hz  -  60  kHz  range  measured  by  the  Noise 
Survey  System. 


1.1  ELF/VLF/LF  Noise 


ELF  /VLF /LF  radio  noise  is  comprised  of  both  man-made  and  natural  signals.  Ex¬ 
amples  of  man-made  signals  are  power  line  harmonics,  communication  signals  and 
interference  from  electrically  powered  machinery;  naturally  occurring  noise  includes 
sferics,  whistlers,  polar  chorus  and  auroral  hiss.  Naturally  occurring  noise  is  dis¬ 
cussed  below,  but  for  a  more  thorough  treatment  see  Helliwell  [24]. 

Sferics  are  typically  the  dominant  source  of  naturally  occurring  low-frequency 
radio  noise.  Even  though  lightning  activity  mainly  occurs  in  tropical  regions  (and 
thus  at  lower  latitudes),  sferics  can  propagate  for  thousands  of  miles  with  little 
attenuation,  so  they  are  seen  in  noise  data  worldwide.  The  amount  of  sferic  activity 
in  a  given  noise  sample  depends  on  the  worldwide  source  distribution  of  lightning 
relative  to  the  receiver  location,  with  nearby  storms  contributing  a  great  deal  and 
distant  storms  contributing  less. 

Whistlers  are  bursts  of  electromagnetic  energy  that  travel  into  the  ionosphere, 
follow  the  lines  of  force  of  the  earth’s  magnetic  field,  and  come  back  into  the  atmo¬ 
sphere  at  roughly  the  same  latitude  and  longitude  in  the  opposite  hemisphere.  The 
dispersion  of  the  medium  in  which  they  travel  causes  them  to  have  a  time-frequency 
characteristic  that  sounds  like  a  “chirp”  or  “whistle”  when  converted  to  an  audible 
signal. 

Polar  chorus  and  auroral  hiss  are  VLF  emissions  that  are  generated  in  the  iono¬ 
sphere  and/or  magnetosphere  and  have  distinct  spectral  signatures.  Both  are  found 
primarily  in  the  polar  regions.  Polar  chorus  refers  to  sequences  of  discrete  signals 
that  often  come  in  bursts  and  have  a  chirping  frequency  characteristic;  auroral  hiss 
is  thermal-like  noise  within  a  given  bandpass  frequency  range.  At  high-latitude  sites, 
it  is  possible  for  auroral  hiss  to  dominate  the  noise  completely  and  drown  out  even 
the  strongest  sferics. 
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Figure  1.1:  Stanford  ELF  Spectrogram,  25  August  1987,  at  00:04  UT. 


It  must  be  noted  that  none  of  the  radio  noise  analyzed  in  this  dissertation  is 
comprised  of  whistlers,  chorus  or  hiss  to  a  noticeable  degree,  so  the  results  contained 
within  are  only  known  to  apply  to  atmospheric  noise. 


1.1.1  Spectral  Content 

A  sample  ELF  spectrogram  from  Stanford,  California,  is  shown  in  Figure  1.1.  The 
first  eight  seconds  include  a  calibration  tone:  harmonics  spaced  10  Hz  apart  with  10 
pico-Tesla  (pT)  magnitudes.  It  is  easily  seen  from  the  calibration  tone  that  notch 
filters  are  installed  at  60,  120,  180  and  300  Hz,  in  order  to  help  remove  power  line 
harmonics  that  otherwise  would  use  up  most  of  the  dynamic  range.  The  rolloff 
above  400  Hz  is  due  to  an  anti-aliasing  filter.  The  vertical  lines  are  sferics,  each 
representing  a  lightning  stroke  that  could  be  many  thousands  of  kilometers  away. 

A  sample  VLF  spectrogram  for  eight  seconds  of  Arrival  Heights  data  is  shown  in 
Figure  1.2.  The  first  one  second  includes  a  calibration  tone:  harmonics  spaced  250 
Hz  apart  with  0.1  pT  magnitudes.  There  is  a  single  pole  highpass  filter  with  a  cutoff 
frequency  of  5  kHz  to  remove  power  line  harmonics.  Once  again  a  number  of  sferics 
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Figure  1.2:  Arrival  Heights  VLF  Spectrogram,  08  May  1995,  at  16:05  UT. 

can  be  seen,  the  strongest  of  which  occur  at  0.7,  6.5  and  7.4  seconds.  The  line  at  10 
kHz  is  an  instrumentation  signal,  and  the  dashed  horizontal  lines  between  10  kHz 
and  14  kHz  on  the  order  of  one  second  in  length  are  from  the  Omega  navigation 
system  [27].  The  horizontal  lines  above  15  kHz  are  from  communication  systems, 
primarily  phase  shift  keyed  digital  systems  transmitting  on  the  order  of  100  -  200 
baud.  When  these  lines  are  crossed  by  the  dark  vertical  lines,  sferics  are  likely  to 
be  causing  bit  errors. 

A  sample  VLF  spectrogram  for  eight  seconds  of  Grafton,  New  Hampshire  data 
is  shown  in  Figure  1.3.  A  great  deal  of  sferic  activity  is  seen,  and  the  sferics  are 
strong  enough  to  clearly  dominate  the  calibration  tone.  This  is  because  Grafton  is 
much  closer  to  thunderstorm  activity  than  Arrival  Heights,  and  in  addition  July  is 
the  peak  of  the  North  American  storm  season  and  00  UT  is  the  peak  of  the  North 
American  diurnal  cycle  [3,  5].  The  term  heavy  sferic  activity  is  used  throughout  this 
dissertation  to  refer  to  the  condition  of  numerous  (and  often  overlapping)  sferics 
as  seen  in  Figure  1.3,  while  the  term  light  sferic  activity  is  used  to  describe  the 
condition  of  sporadic  (but  possibly  large)  impulses  as  in  Figure  1.2. 
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Figure  1.3:  Grafton,  New  Hampshire  VLF  Spectrogram,  20  July  1988,  at  00:05  UT. 

It  is  seen  in  Figures  1.1  to  1.3  that  power  decreases  with  increasing  frequency 
above  a  given  range;  this  is  due  to  the  physical  mechanism  of  lightning  and  to 
propagation  effects.  The  power  spectral  density  of  low-frequency  radio  noise  falls  off 
in  general  as  /“ 2  up  to  2  -  3  kHz,  rises  up  to  10  kHz,  then  drops  off  steeply  (greater 
than  f~4)  throughout  the  VLF  and  LF  ranges  [25,  28].  Propagation  mode  changes 
due  to  the  earth-ionosphere  waveguide  cutoff  account  for  the  2-10  kHz  amplitude 
change. 


1.2  The  Stanford  ELF/VLF  Radio  Noise  Survey 

During  the  years  1985-1986,  eight  ELF/VLF  (10  Hz  -  32  kHz)  radio  noise  mea¬ 
surement  systems,  or  radiometers,  were  installed  at  a  variety  of  high-latitude  and 
mid-latitude  sites  in  an  effort  to  fill  large  gaps  in  the  information  available  on  radio 
noise  in  this  frequency  range  [13,  14].  Other  ELF/VLF  measurement  systems  have 
been  implemented  in  the  past,  but  this  is  the  only  system  of  its  kind  in  terms  of  its 
geographic  coverage  and  continuity  of  simultaneous  data  collection.  The  project  was 
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a  follow-on  to  a  number  of  previous  atmospheric,  ionospheric  and  magnetospheric 
studies  conducted  by  the  Space,  Telecommunications  and  Radioscience  (STAR)  lab¬ 
oratory  at  Stanford  ( e.g Bell  and  Helliwell  [1,  24]). 

The  radiometers  are  located  at  Arrival  Heights,  Antarctica  (AH;  78 °S,  167 °E)-, 
Dunedin,  New  Zealand  (DU;  46°5',  170°#);  Grafton,  New  Hampshire  (GN;  44° N, 
72°W);  Kochi,  Japan  (KO;  33 °N,  1330E);  L’Aquila,  Italy  (AQ;  42° N,  13° E); 
S0ndrestromfj0rd,  Greenland  (SS;  67° N,  51°W);  Stanford,  California  (SU;  37 °N, 
122 °W)\  and  Thule,  Greenland  (TH;  77°  N,  69 °W).  Most  of  the  stations  operated 
much  longer  than  program  expectations,  and  the  systems  at  Stanford  and  Arrival 
Heights  are  still  operating.  They  should  be  able  to  collect  data  beyond  one  solar  cy¬ 
cle.  A  complete  technical  description  of  the  radiometers  is  provided  in  Fraser-Smith 
and  Helliwell  [13],  so  only  an  overview  is  presented  here. 

Each  radiometer  contains  two  receivers,  one  for  the  10  -  400  Hz  frequency  range 
(designated  ELF)  and  the  other  for  the  400  Hz  -  32  kHz  frequency  range  (designated 
VLF).  Each  receiver  has  its  own  pair  of  crossed  loop  antennas,  one  oriented  in  the 
N-S  geomagnetic  direction  and  the  other  in  the  E-W  geomagnetic  direction.  The 
ELF  antennas  are  1164  turn  coils  that  are  either  buried  or  enclosed  in  order  to 
prevent  noise  due  to  wind  induced  motion  of  the  coils  in  the  earth’s  magnetic  field. 
The  VLF  antennas  are  single-turn  triangular  above-ground  loops  18  meters  wide  and 
9  meters  high.  The  Thule,  Greenland  station  has  an  additional  receiver  containing  a 
downconverter  with  a  30  kHz  translation  frequency,  thus  allowing  data  in  the  lower 
part  of  the  LF  range  (30  kHz  -  60  kHz)  to  be  collected  and  analyzed. 

1.2.1  Noise  Amplitude  Data 

Continuous  data  collection  is  obtained  from  the  radiometers  by  monitoring  the  out¬ 
puts  of  a  bank  of  sixteen  narrowband  channel  filters  with  center  frequencies  dis¬ 
tributed  roughly  logarithmically  across  the  10  Hz  -  32  kHz  band.  Each  of  the  32 
filters  (the  N-S  and  E-W  loops  must  be  filtered  separately)  is  a  six  pole  Chebychev 
bandpass  filter  with  a  two-sided  bandwidth  equal  to  five  percent  of  the  center  fre¬ 
quency.  The  sixteen  center  frequencies  and  bandwidths  are  listed  in  Table  1.1:  the 
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Channel 

Frequency 

Bandwidth 

1 

10  Hz 

0.5  Hz 

2 

30 

1.5 

3 

80 

4 

4 

135 

6.75 

5 

275 

13.75 

6 

380 

19 

7 

500 

25 

8 

750  Hz 

37.5 

9 

1  kHz 

50 

10 

1.5 

75 

11 

2 

100 

12 

3 

150 

13 

4 

200 

14 

8 

400 

15 

10.2 

510 

16 

32  kHz 

1600  Hz 

Table  1.1:  Center  frequencies  and  bandwidths  for  the  16  narrowband  channels  of 
the  ELF/VLF  radiometer. 

first  six  are  within  the  ELF  receiver’s  frequency  range  and  the  last  ten  are  within  the 
VLF  receiver’s.  The  10.2  kHz  center  frequency  of  channel  15  was  chosen  specifically 
to  record  signals  from  the  Omega  Navigation  System  (for  propagation  studies),  but 
the  other  center  frequencies  were  chosen  to  be  free  of  man-made  interference. 

Each  filter’s  output  is  passed  through  an  analog  root-mean-square  (RMS)  detec¬ 
tor  that  squares  the  input,  performs  a  time  average,  and  outputs  the  square  root 
of  the  average.  The  RMS  detector  output  is  then  sampled  at  a  rate  of  ten  times 
per  second  by  an  analog  to  digital  converter  and  sent  to  a  digital  computer,  which 
then  computes  the  root-sum-square  of  the  N-S  and  E-W  detector  outputs  to  de¬ 
termine  the  RMS  amplitude  of  the  horizontal  component  of  magnetic  field  for  each 
channel.  The  analog  to  digital  converters  have  a  useful  dynamic  range  of  70  dB,  but 
switchable  gain  amplifiers  in  the  analog  receiver  circuitry  increase  the  total  system 
dynamic  range  to  100  dB.  The  measurements  taken  are  of  absolute  noise  levels. 

In  order  to  save  digital  tape  space,  the  computer  writes  out  only  every  tenth 
sample.  However,  it  also  stores  the  average  and  RMS  values  for  each  minute  (600 
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samples),  along  with  the  minimum  and  maximum  of  the  600  values  for  that  minute. 
The  seasonal  and  diurnal  variations  reported  in  Chapter  2  are  derived  from  these 
one  minute  average  amplitudes. 


1.2.2  Time-Series  Data 

ELF  broadband  time-series  data  are  collected  for  one  minute  every  half  hour,  and 
are  sampled  at  a  1  kHz  sampling  rate  by  a  14  bit  digital  to  analog  converter.  A 
400  Hz  lowpass  filter  is  inserted  before  the  D/A  to  prevent  aliasing.  The  VLF 
and  LF  broadband  time-series  data  are  collected  for  one  minute  every  hour;  the 
sampling  rate  for  these  data  is  62.5  kHz  and  the  D/A  converter  is  16  bits.  The 
anti-aliasing  filter  is  set  at  either  20  kHz  or  30  kHz  depending  on  which  part  of  the 
frequency  spectrum  is  to  be  analyzed.  (The  30  kHz  filter  allows  some  aliasing  at 
lower  frequencies  but  does  not  distort  the  20  kHz  -  30  kHz  range).  The  LF  data  at 
Thule,  once  downconverted,  are  processed  the  same  as  the  VLF  data. 

The  ELF  time-series  data,  like  the  noise  amplitude  data,  are  calibrated  to  abso¬ 
lute  levels  and  digitally  recorded.  The  VLF  time-series  data,  however,  are  recorded 
to  analog  tape  and  digitized  later.  Due  to  concerns  about  the  flatness  of  frequency 
response  in  the  measurement  process,  the  VLF  time-series  data  are  restricted  to  nar¬ 
rowband  analyses  and  their  amplitudes  within  a  given  frequency  band  of  analysis 
are  normalized  to  an  expected  value  of  one. 

Only  one  axis  (the  N-S  loop)  is  recorded  for  the  time  series  data,  so  the  data 
are  collected  with  a  direction  dependent  gain  pattern.  Fortunately,  this  is  shown 
not  to  affect  the  theoretical  results  of  the  dissertation  or  the  comparison  of  these 
results  to  the  noise,  but  it  does  preclude  the  use  of  direction  finding  information  in 
the  statistical  analysis  of  the  noise. 
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1.3  Outline  of  the  Dissertation 

This  dissertation  is  concerned  with  the  statistical  analysis  and  modeling  of  low- 
frequency  radio  noise  and  the  improvement  of  low-frequency  communications.  Chap¬ 
ter  2  discusses  statistical  analysis  results,  Chapters  3  and  4  are  concerned  with  noise 
modeling,  and  Chapter  5  addresses  low-frequency  communications. 

The  statistics  analyzed  in  Chapter  2  include  seasonal  and  diurnal  variations, 
amplitude  probability  distributions  (APDs),  impulse  interarrival  time  distributions, 
correlations  between  noise  impulses,  and  background  noise  statistics.  Noise  char¬ 
acteristics  derived  from  the  analyses  are  presented,  and  are  used  in  the  following 
chapters  to  aid  in  the  modeling  problem. 

In  Chapter  3,  APDs  derived  from  thousands  of  hours  of  ELF/VLF/LF  noise 
survey  data  are  compared  to  the  APDs  of  various  low-frequency  noise  models.  The 
most  accurate  model  for  representing  low-frequency  noise  APDs  is  found  to  depend 
on  geographic  location,  time  of  year  and  day,  bandwidth,  and  center  frequency, 
but  two  of  the  simplest  models  ( i.e .,  each  with  only  two  parameters)  are  found  to 
give  extremely  good  performance  in  general.  These  are  the  Hall  and  alpha-stable 
(or  c*-stable)  models,  both  of  which  approximate  the  Rayleigh  distribution  for  low 
amplitude  values  but  decay  with  an  inverse  power  law  for  high  amplitude  values. 
It  is  concluded  that  the  Hall  model  is  the  optimal  choice  in  terms  of  accuracy  and 
simplicity  for  locations  exposed  to  heavy  sferic  activity  ( e.g .,  low-latitude  regions), 
and  the  o-stable  model  is  best  for  locations  relatively  distant  from  heavy  sferic 
activity  (e.g.,  high-latitude  regions). 

Chapter  4  develops  a  new  clustering  Poisson  atmospheric  noise  model  based  on 
the  statistical  characteristics  of  the  noise  data  found  in  Chapters  2  and  3.  This 
model  is  similar  to  several  previously  known  statistical-physical  models,  but  in  ad¬ 
dition  accounts  for  the  clustering  of  sferics.  It  is  verified  that  the  clustering  model 
accurately  predicts  the  statistical  features  found  in  low-frequency  noise  data. 

Chapter  5  compares  the  bit  error  rate  (BER)  performance  of  receivers  specifically 
designed  for  a-stable  noise  with  that  of  conventional  nonlinear  analog  receivers  used 
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with  impulsive  noise,  such  as  clippers,  limiters  and  log-correlators.  The  communi¬ 
cation  signal  formats  examined  are  quadrature  phase-shift  keying  (QPSK)  and  16 
point  quadrature  amplitude  modulation  (16  QAM).  Hundreds  of  simulations  using 
time-series  data  from  various  times  and  locations  and  at  various  center  frequencies 
and  bandwidths  are  performed,  and  the  following  results  are  found  uniformly:  for 
QPSK  signals,  virtually  no  performance  improvement  is  gained  when  using  an  te¬ 
stable  receiver  instead  of  the  best  conventional  receiver,  but  for  16  QAM  signals,  an 
improvement  of  several  dB  is  gained  by  using  an  a-stable  receiver. 

Finally,  Chapter  6  summarizes  the  main  results  of  this  dissertation  and  suggests 
topics  of  future  research. 


1.4  Contributions 

Here  is  a  brief  summary  of  the  contributions  made  in  this  dissertation. 

•  Statistical  analysis  of  thousands  of  hours  of  globally  collected  low-frequency 
radio  noise  (Chapter  2). 

•  Comparison  of  APD  models  used  for  modeling  low-frequency  noise  (Chap¬ 
ter  3). 

•  Development  of  a  clustering  Poisson  model  for  atmospheric  noise  and  proof 
that  the  model  accurately  predicts  low-frequency  noise  statistics  (Chapter  4). 

•  Comparison  of  receivers  designed  for  a-stable  noise  to  conventional  receivers 
(Chapter  5). 


Chapter  2 


Measured  Statistics 


2.1  Introduction 

This  chapter  presents  an  overview  of  some  of  the  statistical  characteristics  of  nat¬ 
urally  occurring  low-frequency  radio  noise.  The  statistical  properties  discussed 
include  long-term  noise  averages,  amplitude  probability  distributions,  impulse  in¬ 
terarrival  time  distributions,  correlations  between  impulses,  and  background  noise 
statistics. 

2.2  Long-Term  Noise  Averages 

The  noise  amplitude  data  described  in  Section  1.2.1  are  collected  continuously  at 
the  radiometer  sites  over  long  periods  of  time,  and  their  absolute  noise  levels  are 
well  calibrated.  It  is  thus  possible  to  derive  statistically  significant  seasonal  and 
diurnal  variations  of  ELF /  VLF  noise  from  them.  Since  the  noise  is  primarily  caused 
by  lightning  occurring  throughout  the  world,  these  variations  are  related  to  seasonal 
weather  patterns  and  global  climate  change  [15].  This  section  provides  an  overview 
of  the  seasonal  and  diurnal  variation  calculations  performed  with  the  data;  more 
extensive  results  are  contained  in  [3,  4,  5]. 

To  provide  a  basic  context  for  the  ELF/VLF  noise  amplitude  measurements, 
Figure  2.1  shows  the  one  hour  average  noise  amplitudes  over  the  course  of  one 
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Arrival  Heights,  Antarctica  JUN  94  Noise  Survey  Average 


Figure  2.1:  Noise  averages  in  the  10  Hz  frequency  band  recorded  June  1994,  at 
Arrival  Heights,  Antarctica.  Each  of  the  720  points  on  this  plot  is  an  average  of  one 
hour  of  data. 

month  for  one  channel,  the  10  Hz  band  measured  at  Arrival  Heights  during  the 
month  of  June  1994.  Each  of  the  720  points  on  this  graph,  representing  one  hour, 
is  an  average  of  roughly  32,000  noise  filter  output  samples  (not  36,000  because  of 
calibration  periods).  The  data  consist  of  both  random  and  diurnal  variations;  they 
sometimes  show  occasional  short  duration  impulses  due  to  both  man-made  and 
natural  interference  as  well.  The  entire  database  contains  thousands  of  these  plots, 
one  for  each  station,  month  and  channel. 

The  unit  fT/\/Hz  is  essentially  the  square  root  of  power  spectral  density,  ob¬ 
tained  in  this  case  by  dividing  the  RMS  filter  amplitude  output  (in  fT)  by  0.707, 
the  square  root  of  the  0.5  Hz  bandwidth  for  the  10  Hz  channel  filter.  The  data  are 
presented  in  fT  because  the  system  detects  magnetic  field.  The  vertical  electric  field 
component  and/or  the  power  of  the  incoming  signal  may  be  obtained  using  377  0 
as  the  impedance  of  free  space,  but  this  is  an  approximation  (albeit  usually  a  good 
one)  that  assumes  the  impinging  electromagnetic  waves  are  planar.  To  convert  to 
electric  field  under  this  assumption,  the  relation  B  =  y/Ji^E  =  E/c  may  be  used  to 
determine  that  1  fT  is  equivalent  to  0.300  ^V/m.  If  it  is  desired  to  relate  magnetic 
field  (or  magnetic  flux  density)  B  to  magnetic  intensity  H,  the  relation  B  =  fj,0H 
can  be  used  to  find  that  1  fT  is  equivalent  to  7.958  x  10-4  ^A/m. 

2.2.1  Seasonal  Variations 

Determining  the  seasonal  variations  of  natural  radio  noise  over  the  course  of  many 
years  requires  two  steps:  (i)  monthly  averages  are  computed  for  each  station,  month, 
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year,  and  frequency  band  (i.e.,  all  the  data  in  Figure  2.1  are  averaged  to  come  up 
with  one  number  for  June  1994  at  Arrival  Heights  in  the  10  Hz  frequency  band), 
and  (ii)  all  the  years  are  then  averaged  together,  resulting  in  an  overall  variation  by 
month  at  each  station  and  in  each  band. 

Figures  2.2  and  2.3  show  the  monthly  variations  of  the  noise  level  for  all  channels 
at  Stanford,  California,  and  the  error  bars  show  the  standard  deviations  from  year 
to  year,  i.e.,  a  small  error  bar  indicates  little  variation  from  one  year  to  the  next  for 
that  month.  Note  that  the  error  bars  are  largely  unrelated  to  the  standard  deviation 
of  individual  noise  samples,  which  can  be  quite  large;  they  are  also  unrelated  to  the 
standard  deviation  of  the  total  average,  which  is  minute  since  each  point  on  these 
plots  is  an  average  of  millions  of  sample  values.  Note  that  each  individual  graph  has 
its  own  scale.  All  24  hours  of  each  day  are  included;  for  seasonal  variations  using 
only  specific  diurnal  time  periods,  refer  to  [3]. 

A  strong  annual  variation  with  a  peak  in  the  northern  hemisphere  summer  can 
be  seen  in  most  of  the  frequency  channels  at  Stanford.  Although  not  shown,  it 
exists  at  Spndrestrpm  as  well.  Arrival  Heights  and  Dunedin,  however,  have  some 
variations  with  a  peak  in  the  northern  hemisphere  winter.  The  differences  between 
annual  variations  among  stations  and  from  channel  to  channel  can  be  attributed  to 
different  patterns  in  local  and  distant  lightning,  with  the  higher  frequency  variations 
being  influenced  more  strongly  by  closer  sources.  The  seasonal  variation  of  global 
lightning  is  such  that  southern  hemisphere  locations  are  generally  more  active  in 
the  northern  hemisphere  winter  and  northern  hemisphere  locations  are  more  active 
in  the  northern  hemisphere  summer  [4,  19]. 

In  Figure  2.3,  for  frequencies  in  the  range  1  -  1.5  kHz,  the  error  bars  are  too 
large  relative  to  the  data  to  extract  a  statistically  significant  seasonal  variation.  This 
band  is  within  the  range  of  the  earth-ionosphere  waveguide  cutoff  frequencies,  where 
noise  does  not  propagate  far  and  the  receiver  predominantly  picks  up  the  fields  from 
local  sources.  In  addition,  the  32  kHz  band  at  Stanford  is  known  to  be  contaminated 
by  man-made  noise,  resulting  in  a  weak  variation  with  large  error  bars. 

At  the  other  frequencies,  however,  a  monthly  variation  strongly  dominates  the 
uncertainty  of  the  error  bars.  Given  the  number  of  data  points  and  years  involved, 
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Figure  2.2:  Monthly  variation  of  ELF/VLF  radio  noise  at  Stanford,  California,  for 
the  eight  lowest-frequency  channels.  The  years  1986  to  1993  are  included. 
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Figure  2.3:  Monthly  variation  of  ELF/VLF  radio  noise  at  Stanford,  California,  for 
the  eight  highest-frequency  channels.  The  years  1986  to  1993  are  included. 
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this  reflects  a  definite  seasonal  trend  in  source  quantity  ( i.e .  lightning  flash  rate) 
and/or  location.  It  is  shown  in  [4]  that  the  seasonal  variations  seen  in  the  data 
correlate  strongly  with  estimates  of  global  lightning  flash  rates. 

2.2.2  Diurnal  Variations 

Figures  2.4  and  2.5  show  diurnal  variations  of  the  noise  level  for  January  and  July  at 
Arrival  Heights,  Antarctica,  for  the  eight  lowest  frequency  channels.  The  error  bars 
indicate  the  standard  deviations  from  year  to  year,  i.e.,  small  error  bars  indicate  little 
variation  from  one  year  to  the  next  in  the  diurnal  cycle  of  a  particular  month.  Once 
again,  the  error  bars  are  largely  unrelated  to  the  standard  deviation  of  individual 
noise  samples  or  to  the  standard  deviation  of  the  total  average. 

Changes  in  average  noise  levels  from  year  to  year  for  each  month  are  removed 
from  the  error  bar  calculations,  since  otherwise  they  would  artificially  enlarge  the 
diurnal  variation  error  bars.  The  normalization  is  performed  in  three  parts:  (i)  for 
each  month  and  year,  a  noise  average  for  the  entire  month  is  computed  to  produce  a 
single  value,  (ii)  for  each  month,  the  resulting  values  from  (i)  for  the  different  years 
are  averaged  together  to  give  one  total  average  reference  value  per  month,  and  (iii) 
all  the  data  are  normalized  (by  subtracting  the  difference  between  the  corresponding 
values  from  parts  (i)  and  (ii))  such  that  differences  in  total  monthly  averages  from 
year  to  year  are  removed,  i.e.,  each  total  monthly  average  now  equals  that  month’s 
reference  value.  Thus  the  error  bars  are  truly  an  indication  of  the  variation  of  the 
diurnal  cycle.  In  most  cases  the  error  bars  are  small  compared  to  their  respective 
data,  indicating  little  variation  of  diurnal  cycles  from  year  to  year. 

Arrival  Heights  (at  77.8°  S  latitude)  sees  roughly  equal  contributions  from  storms 
at  all  longitudes,  especially  at  the  lowest  frequencies.  The  10  Hz  band  thus  exhibits 
a  diurnal  pattern  in  phase  with  the  overall  worldwide  distribution,  a  broad  peak 
from  14-22  UT,  as  clearly  seen  in  Figure  2.4.  From  June  to  September,  however, 
the  peak  shifts  to  roughly  22-00  UT,  as  seen  in  Figure  2.5.  This  indicates  a  strong 
contribution  from  North  American  storms  during  their  peak  season  [15]. 

Although  not  shown,  large  diurnal  variations  are  seen  in  most  of  the  frequency 
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Figure  2.4:  Diurnal  variation  of  ELF /VLF  radio  noise  at  Arrival  Heights,  Antarctica, 
during  the  month  of  January  for  the  eight  lowest-frequency  channels.  The  years  1985 
to  1994  are  included. 
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Arrival  Heights,  Antarctica,  JUL  Diurnal  Variation  (fT/VHz) 
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Figure  2.5:  Diurnal  variation  of  ELF /VLF  radio  noise  at  Arrival  Heights,  Antarctica, 
during  the  month  of  July  for  the  eight  lowest-frequency  channels.  The  years  1985 
to  1994  are  included. 
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channels  at  most  of  the  stations  [5].  However,  the  phases  of  these  diurnal  variations 
can  depend  strongly  on  month,  frequency,  and  especially  station.  Variation  by 
station  is  primarily  due  to  the  diurnal  signature  of  global  lightning,  for  which  it  is 
known  that  lightning  over  North  America  peaks  at  roughly  00  UT,  lightning  over 
South  America  peaks  at  roughly  20  UT,  lightning  over  Europe  and  Africa  peaks  at 
roughly  16  UT,  and  lightning  over  Southeast  Asia  peaks  at  roughly  08  UT  [15,  19]. 
The  differences  in  diurnal  variations  with  respect  to  frequency  can  be  attributed  to 
different  patterns  in  local  and  distant  lightning,  with  the  higher  frequency  variations 
being  influenced  more  strongly  by  closer  sources. 


2.3  Amplitude  Probability  Distributions 

The  most  measured  and  modeled  statistic  of  low-frequency  radio  noise,  next  to 
absolute  power  level,  is  the  first-order  APD.  This  section  gives  an  overview  of  APD 
characteristics;  APD  modeling  is  covered  in  Chapters  3  and  4. 

Since  wideband  time-series  data  generally  contain  significant  man-made  inter¬ 
ference,  the  noise  APD  is  usually  analyzed  within  a  relatively  narrow  frequency 
range.  In  such  narrowband  analysis,  the  broadband  time-series  data  are  digitally 
downconverted  (i.e.,  frequency  translated)  using  various  center  frequencies,  and  low- 
pass  filtered  using  various  cutoff  frequencies.  Statistics  are  then  derived  from  the 
resulting  lowpass  equivalent  signals. 

A  narrowband  signal  n(t )  can  be  written 

n{t)  =  n/(f )  cos(2tt ft)  —  riQ(t )  sin(27r/i),  (2.1) 

where  nj(t)  is  the  inphase  component,  n<p(t)  is  the  quadrature  component,  and  / 
is  the  center  frequency.  The  signals  nj(f)  and  nQ(t)  are  then  lowpass  signals  with 
bandwidth  much  smaller  than  /.  The  complex-analytic  representation  is  written 


n(f)  =  A{t)ejV”ft+&{t)\ 


(2.2) 
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where  the  magnitude  (or  envelope)  A(t)  is 

A(t)  =  +  nl(t),  (2.3) 

and  the  phase  0(t)  is 

0(t)  =  arctan  (2.4) 

Using  this  latter  notation  the  signal  n(t)  can  also  be  written 

n(t)  =  R[A(t)ej(2*ft+&(%  (2.5) 

where  3?()  indicates  the  real  part  of  the  argument. 

The  noise  envelope  A(t)  is  a  random  process,  and  its  first-order  statistics  at  a 
given  time  are  specified  by  its  APD.  The  APD  is  defined  as  the  probability  that  the 
noise  envelope  A  takes  on  a  value  larger  than  some  given  value  a:  P(A  >  a).  Other 
commonly  used  statistical  definitions  that  characterize  A(t)  are  (i)  the  cumulative 
distribution  function  (CDF),  FA(a ),  which  is  one  minus  the  APD,  (ii)  the  probability 
density  function  (pdf),  /4(a),  which  is  the  derivative  of  the  CDF,  and  (iii)  the 
voltage  deviation  V^,  defined  as  10  log (E[A?\/ E2\A\),  which  serves  as  an  indicator  of 
the  impulsiveness  of  the  noise.  The  phase  Q(t)  of  atmospheric  noise  has  long  been 
known  to  have  a  distribution  that  is  uniform  over  the  angles  -7 r  to  tt  [33];  several 
checks  on  the  noise  survey  data  provide  additional  evidence  that  this  is  true. 

Since  low-frequency  noise  is  thought  of  as  Gaussian  noise  plus  an  impulsive  noise 
component,  its  pdf  looks  like  a  bell  curve  but  with  heavier  tails.  Likewise,  since 
the  envelope  distribution  of  Gaussian  noise  is  Rayleigh,  the  envelope  distribution 
of  low-frequency  noise  appears  Rayleigh  but  with  a  heavier  tail.  The  difference 
between  a  typical  low-frequency  noise  envelope  pdf  and  the  Rayleigh  pdf  is  shown 
in  Figure  2.6,  for  July  1986  data  from  Thule,  Greenland,  in  the  35.6  -  37.6  kHz 
range.  Noise  envelope  pdf’s  are  found  to  decay  with  an  inverse  power  law  for  large 
values  (out  to  some  limit  set  by  the  dynamic  range  of  the  system),  so  the  Rayleigh 
pdf  has  a  tail  that  rolls  off  too  quickly  to  accurately  represent  the  noise  data  pdf. 
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Figure  2.6:  Thule  envelope  pdf  compared  to  a  Rayleigh  distribution  with  the  same 
initial  slope.  The  data  are  normalized  such  that  E[A\  =  1. 

2.4  Impulse  Interarrival  Distributions 

Thousands  of  hours  of  ELF/VLF/LF  data  were  analyzed  to  determine  impulse 
interarrival  distributions.  In  each  case,  a  set  of  threshold  values  (and  their  negatives) 
were  chosen  such  as  in  Figure  2.7,  a  plot  of  52  seconds  of  ELF  broadband  time-series 
data.  For  each  occurrence  of  an  upward  crossing  of  a  given  positive  threshold  value 
or  a  downward  crossing  of  the  respective  negative  threshold  value,  the  following  are 
noted:  (i)  time  of  occurrence  of  the  threshold  crossing,  linearly  extrapolated  from 
the  two  samples  between  which  the  level  is  crossed,  (ii)  the  maximum  (or  minimum) 
extent  of  the  sferic,  and  (iii)  the  amount  of  time  the  noise  remains  above  (or  below) 
the  threshold  crossing  ( i.e the  duration,  also  linearly  extrapolated  from  the  two 
closest  sample  values).  Crossings  less  than  one  millisecond  apart  are  assumed  to 
be  from  the  same  impulsive  waveform  and  are  ignored.  For  a  given  +/—  threshold 
value,  all  positive  and  negative  sferics  are  treated  together,  and  the  time  differences 
are  determined  from  one  sferic  to  the  next.  The  result  is  a  set  of  impulse  spacings, 
or  interarrival  times. 
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Figure  2.7:  Stanford  ELF  broadband  data,  01  April  1990,  at  01:04  UT.  Dotted  lines 
are  shown  at  ±5  and  ±10  times  the  norm  of  the  noise. 
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Figure  2.8:  Impulse  interarrival  pdf  for  Arrival  Heights,  December  1995,  for  the 
25.5  -  27.5  kHz  noise  band,  with  the  impulse  detection  threshold  set  at  20  times  the 
norm  of  the  noise. 
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The  set  of  successive  interarrival  times  Tj(fc)  form  a  discrete  random  process. 
The  marginal  pdf  of  Tj(k)  is  assumed  not  to  depend  on  the  index  k,  and  so  can  be 
found  from  a  histogram  of  the  interarrival  times  of  the  data.  For  every  data  sample 
analyzed,  the  resulting  interarrival  time  pdf  fj ,(t)  was  found  to  be  of  the  form  of 
Figure  2.8;  this  figure  shows  /x/( t)  for  Arrival  Heights  December  1995  data  in  the 
25.5  -  27.5  kHz  noise  band,  with  the  impulse  threshold  set  at  20  times  the  norm  of 
the  noise. 

Since  Figure  2.8  has  a  logarithmic  ordinate  axis,  functions  of  the  form  ce~Xt 
appear  as  straight  lines  with  slope  -A  and  intercept  log  c.  Thus  an  exponential  pdf, 
indicative  of  a  simple  Poisson  process,  would  appear  as  a  straight  line.  In  all  the 
data  and  for  all  threshold  values,  an  approximately  straight  line  is  found  for  times 
greater  than  roughly  one  second.  For  values  of  Tj  between  zero  and  one,  however,  the 
slope  is  more  steeply  downward  curving.  This  implies  clustering,  since  it  indicates 
that  many  of  the  impulses  are  bunched  together  with  short  times  between  them. 
Section  4.4  shows  that  a  clustering  Poisson  assumption  predicts  the  behavior  of  the 
pdf  fait)  very  well. 


2.5  Correlations  Between  Impulses 


Since  the  magnitude  and  direction  of  each  sferic’s  threshold  crossing  is  noted,  it  is 
possible  to  check  for  correlations  between  maximum  amplitude  levels  of  adjacent 
impulses  and  also  between  maximum  amplitudes  and  interarrival  times.  In  all  the 
data  samples  examined  and  at  a  variety  of  threshold  levels,  however,  no  significant 
correlations  were  found  in  adjacent  amplitude  levels  or  between  amplitudes  and 
interarrival  times,  validating  certain  independence  assumptions  of  the  model  to  be 
proposed  in  Section  4.2.2.  The  correlation  coefficient  between  adjacent  interarrival 
times  is  roughly  .15  -  .20. 
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2.6  Background  Noise  Statistics 

As  mentioned  previously,  it  has  long  been  held  that  atmospheric  noise  can  be  viewed 
as  a  superposition  of  Gaussian  noise  and  distinct  impulses  [12,  17].  The  Gaussian 
component  is  the  background  noise,  which  results  from  the  superposition  of  numerous 
low-level  sources,  including  distant  sferics.  The  impulses  are  caused  by  lighting  from 
relatively  nearby  thunderstorms. 

For  narrowband  noise  with  bandwidth  B,  the  background  noise  samples  are  com¬ 
monly  modeled  as  independent  if  they  are  spaced  at  the  Nyquist  rate  1/2B.  This 
assumption  about  the  noise  survey  data  is  tested  by  removing  all  the  noise  impulses 
that  stand  out  above  the  background  noise  level  and  performing  correlation  statistics 
on  the  samples  that  are  left  (usually  more  than  95%  of  them).  It  is  found  that  the 
resulting  background  noise  is  indeed  well  modeled  as  white  Gaussian  noise  (WGN), 
albeit  with  the  tails  chopped  off.  This  characteristic  is  used  in  the  receiver  design 
problem  of  Chapter  5. 


2.7  Conclusions 

This  chapter  sets  the  framework  for  the  rest  of  the  dissertation  by  introducing  some 
of  the  statistical  properties  of  low-frequency  radio  noise.  These  properties  are  used 
in  the  following  chapters  to  derive  both  models  of  the  noise  and  receiver  structures 
for  communicating  in  it.  The  primary  results  of  this  chapter  are: 

•  Seasonal  and  Diurnal  variations  correlate  well  with  global  lightning  patterns 
and  can  be  used  to  track  global  climate  change. 

•  The  narrowband  noise  envelope  pdf  of  low-frequency  noise  has  the  functional 
form  of  a  Rayleigh  pdf  at  low  values  of  dynamic  range  but  decays  with  an 
inverse  power  law  for  high  values. 


•  Impulse  interarrival  time  distributions  indicate  clustering  of  sferics  with  bursts 
on  the  order  of  one  second  in  length. 
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•  No  appreciable  correlations  are  found  between  the  maximum  amplitudes  of 
adjacent  sferics  or  between  these  maximum  amplitudes  and  interarrival  times. 

•  Low-frequency  noise  can  be  viewed  as  the  superposition  of  white  Gaussian 
noise  and  clustered  noise  impulses  from  sferics. 


Chapter  3 

Comparison  of  APD  Models 


3.1  Introduction 


It  is  stated  in  Chapter  2  that  one  of  the  most  commonly  modeled  statistics  in 
atmospheric  noise  studies  is  the  noise  envelope  APD,  or  equivalently  the  pdf.  In 
this  chapter,  the  first-order  amplitude  probability  distributions  of  ELF/VLF/LF 
noise  data  are  determined  and  compared  to  several  noise  models  used  in  recent 
literature:  the  Hall  model  [23],  the  Field- Lewinstein  (or  F-L)  model  [12]  and  the 
o-stable  model  [40].  These  were  chosen  over  a  number  of  other  noise  models  because 
of  their  popularity,  accuracy  and  relative  simplicity. 

For  a  given  data  histogram,  each  model’s  parameters  are  adjusted  to  minimize 
the  expected  value  of  the  log  error  squared  between  the  histogram  and  the  model’s 
estimate  of  it.  This  log  error  metric,  called  the  mean-square  log  error  (MSLE),  is 
defined  as 


MSLE  = 


dx, 


(3.1) 


where  f(x)  is  the  data  pdf  and  f(x)  is  the  model’s  estimate  of  it.  Note  that  the 
expression  (3.1)  is  similar  to  the  relative  entropy  definition  in  information  theory  [6], 
except  for  squaring  the  term  in  parentheses. 

After  optimizing  the  parameters  of  each  model,  the  minimum  errors  achieved 
for  the  individual  models  are  compared  and  the  best  model  is  determined.  The 
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best  model  depends  on  geographic  location,  time  of  year  and  day,  frequency,  and 
bandwidth,  but  the  Hall  and  a-stable  models  are  found  to  give  extremely  good 
performance  in  general.  Both  models  have  Rayleigh  characteristics  for  low  amplitude 
values  but  decay  with  an  inverse  power  law  for  high  amplitude  values.  The  results 
show  that  the  Hall  model  is  the  optimal  choice  in  terms  of  accuracy  and  simplicity 
for  locations  exposed  to  heavy  sferic  activity  (e.g.,  low-latitude  regions),  and  the 
a-stable  model  is  best  for  locations  relatively  distant  from  heavy  sferic  activity 
(e.g.,  high-latitude  regions).  In  addition,  the  commonly  used  voltage  deviation  (Vd) 
parameter  is  found  to  be  only  weakly  related  to  the  amount  of  sferic  activity  in  the 
noise  waveform. 

3.2  Amplitude  pdf  Models 

Extensive  research  has  been  conducted  for  over  forty  years  in  order  to  characterize 
and  model  the  first-order  APDs  of  atmospheric  radio  noise  (e.g.,  [7,  33,  36,  47]). 
Some  APD  models  are  based  entirely  on  intuitive  reasoning  and/or  fitting  the  data 
to  mathematical  functions  (these  are  called  empirical  models);  others  start  with 
assumptions  on  noise  source  distributions  and  the  propagation  of  noise  impulses  to 
the  receiver  (statistical-physical  models).  Empirical  models  in  general  have  been 
based  on  crude  curve  fittings  to  limited  data,  but  nonetheless  a  number  of  good 
models  have  been  developed,  including  the  Hall  and  Field- Lewinstein  models.  They 
are  typically  simpler  and  more  mathematically  tractable  than  other  models,  but 
their  parameters  are  often  unrelated  to  the  physical  processes  that  create  the  noise. 

Statistical-physical  models  take  into  account  the  underlying  physical  processes 
of  impulsive  noise,  but  usually  are  difficult  to  work  with  mathematically.  Making 
them  tractable  often  requires  approximations  that  are  known  not  to  be  true  for 
atmospheric  noise,  such  as  assuming  the  impulsive  sources  are  distributed  indepen¬ 
dently  and  uniformly  in  space  and  time.  Nonetheless,  statistical-physical  models 
prove  useful  in  many  applications. 

The  most  widely  known  statistical-physical  models  are  the  class  A  and  B  noise 
models  developed  by  Middleton  [33].  These  are  not  considered  explicitly  in  this 
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dissertation  due  to  their  complexity,  but  the  class  B  model  proves  to  be  very  accurate 
since  it  is  a  generalization  of  both  the  Hall  and  a-stable  models.  It  should  be  noted 
that  class  A  noise  is  defined  for  cases  in  which  the  bandwidth  of  the  input  noise  is 
comparable  to  or  less  than  the  detection  bandwidth,  and  class  B  noise  is  defined  for 
cases  where  the  noise  bandwidth  is  greater  than  the  bandwidth  of  the  detector  (i.e. 
impulsive  noise  inputs  produce  transients  in  the  receiver).  It  is  the  latter  case  that 
applies  in  this  study,  especially  at  lower  latitudes,  since  the  sferics  in  the  noise  data 
have  much  wider  bandwidths  than  those  used  to  analyze  the  noise. 

The  complete  set  of  noise  model  distributions  from  which  the  Hall,  a-stable 
and  F-L  were  selected  includes  power-Rayleigh  (or  Weibull),  Laplace,  lognormal, 
hyperbolic,  and  any  promising  mixture  processes  or  piecewise  combinations  of  these, 
both  with  and  without  the  Gaussian  or  Rayleigh  distribution  included  as  well.  This 
encompasses  most,  if  not  all,  of  the  commonly  known  APD  models  ( e.g [16, 18,  34]). 


3.2.1  Hall  Model 

The  Hall  Model  is  presented  and  explained  extensively  in  [23],  but  only  the  first 
order  pdf  of  the  narrowband  envelope  A  is  needed  here.  The  Hall  model  specifies 
this  envelope  pdf  /4(a)  as  the  two  parameter  distribution 

fA(a)  =  (m-  l)7m_1  a21m±.,  a  >  0,  (3.2) 

called  the  Hall  pdf  from  this  point  on.  The  term  7  is  a  scaling  factor  and  the  term 
m  determines  the  impulsiveness  of  the  noise.  Note  that  A  has  infinite  variance  for 
m  <  3. 


3.2.2  Field-Lewinstein  Model 

The  Field-Lewinstein  model  is  an  empirical  model  developed  from  the  assumption 
that  atmospheric  noise  is  composed  of  impulsive  noise  superimposed  on  a  back¬ 
ground  of  low-level  Gaussian  noise.  The  envelope  A  is  approximated  as  the  sum  of 
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a  Rayleigh  distributed  random  variable  (for  the  Gaussian  component)  and  a  power- 
Rayleigh  distributed  random  variable  (for  the  impulsive  component),  so  the  envelope 
pdf  is  the  convolution  of  the  two  densities: 


Ia  (a) 


aa 


a- 1 


Ra 


a  >  0. 


(3.3) 


Equation  (3.3)  is  referred  to  as  the  F-L  distribution  from  this  point  on.  Note  that 
it  is  a  three  parameter  distribution  with  no  closed  form,  so  it  is  somewhat  difficult 
to  work  with  mathematically. 


3.2.3  Alpha-Stable  Model 

The  a-stable  pdf  is  directly  specified  in  the  characteristic  function  domain.  The 
characteristic  function  $x(u;)  of  a  random  variable  X  is  essentially  a  Fourier  trans¬ 
form  of  the  pdf: 

=  E[e^x)  =  f"  fx{x)ei"*dx, 

J  —  OO 

and  the  a-stable  characteristic  function  is  defined  as 

*jr(w)  =  e"7M“, 

a  two  parameter  distribution.  The  general  form  of  the  distribution  includes  two 
more  parameters,  defining  an  absolute  shift  and  a  skew,  but  these  may  be  eliminated 
by  assuming  the  noise  is  distributed  symmetrically  about  zero.  (This  the  case  for 
atmospheric  noise  time-series  data). 

For  a  =  2  the  characteristic  function  defines  a  Gaussian  distribution  with  mean 
zero  and  variance  27 ;  for  a  =  1  it  defines  a  Cauchy  distribution  with  parameter  7. 
Thus  the  Cauchy  and  Gaussian  distributions  are  forms  of  the  a-stable  distribution. 
For  a  close  to  two  the  distribution  is  essentially  Gaussian  but  with  heavier  tails. 
The  a-stable  distribution  for  7  =  1  and  various  values  of  a  is  shown  for  positive 
values  in  Figure  3.1.  A  significant  difference  in  tail  rolloff  is  seen  even  for  a  =  1.95 
as  compared  to  the  Gaussian  case,  a  =  2.0. 
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Figure  3.1:  The  a-stable  pdf  for  a  =  1.0  (Cauchy),  1.6,  1.95  and  2.0  (Gaussian), 
and  7  =  1.0;  only  positive  abscissa  values  are  shown. 


If  narrowband  atmospheric  noise  is  a-stable  distributed,  the  inphase  and  quadra¬ 
ture  components  are  jointly  a-stable  with  characteristic  function 

$Xi,xQ(tn)  =  Ele^X’+^x^]  =  e-^'Q,  (3.4) 

where  u  =  (wi,u?2).  The  corresponding  envelope  distribution  is  the  Fourier-Bessel 
transform  of  e-7pQ: 


roo 

/yi(a)  =  a  pe~'yp  J0(ap)dp,a  >  0.  (3.5) 

Jo 

Note  for  a  =  2  that  this  is  the  Rayleigh  distribution.  For  a  close  to  two  it  is 
essentially  a  Rayleigh  distribution  but  with  a  heavier  tail.  From  this  point  on,  the 
/a (a)  of  (3.5)  will  be  referred  to  as  the  a-stable  envelope  pdf. 

The  a-stable  pdf  does  not  exist  in  closed  form  except  for  a  =  1  or  a  =  2,  but  it 
can  be  approximated  numerically  without  a  great  deal  of  computational  complexity. 
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In  addition,  the  parameter  7  is  a  simple  scaling  factor  such  that 

h :(*;«>  7)  = 

a  'Ja 

and  likewise  for  the  envelope  pdf 

/^(o,  01,7)  t-/a( 

ot  ot 

so  a  lookup  table  must  vary  only  over  the  one  parameter  a.  Such  a  lookup  table 
method  is  used  in  these  analyses. 

The  a-stable  model  is  an  empirical  model,  but  like  the  Hall  model  it  does  have 
some  physical  justification.  Nikias  and  Shao  [37]  show  that  under  certain  assump¬ 
tions  on  the  underlying  noise  sources  and  the  propagation  characteristics  between 
them  and  the  receiver,  atmospheric  noise  is  expected  to  exhibit  an  a-stable  pdf. 
These  assumptions,  such  as  defining  sources  to  be  independently  distributed  in  space 
and  time,  are  not  true  in  practice;  however,  they  are  approached  closely  enough  in 
some  cases  to  explain  the  accuracy  of  the  a-stable  model. 


3.3  Determining  APDs  of  the  data 

For  data  analysis,  the  underlying  probability  distribution  /a  (a)  of  the  noise  envelope 
must  be  derived  from  the  data,  so  it  is  assumed  that  the  desired  ensemble  statistics 
can  be  determined  using  corresponding  time  averages  (for  time  intervals  of  interest). 
The  data  consist  of  one  minute  segments  spaced  one  hour  apart,  however,  and  it  is 
known  that  the  noise  is  not  stationary  because  the  mean  and  RMS  values  have  a 
diurnal  and  seasonal  variation.  To  compensate  for  the  shift  in  absolute  amplitude 
levels  between  one  minute  segments,  a  normalization  is  performed  on  each  segment 
such  that  the  norm  of  A  equals  one,  i.e.,  F[|A|]  =  E[A\  =  1. 

A  normalization  in  amplitude  is  used  rather  than  the  more  common  normaliza¬ 
tion  in  power,  E[A2]  =  1,  because  of  the  large  effect  a  few  outliers  (large  impulsive 
values)  can  have  on  the  second  moment.  When  normalizing  to  the  second  moment, 
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just  a  few  large  sferics  in  one  sample  compared  to  another  can  shift  the  distribu¬ 
tions  with  respect  to  each  other  even  though  their  statistics  are  otherwise  similar. 
When  using  normalization  in  amplitude,  however,  data  taken  over  the  course  of 
hours  and/or  days  have  histograms  that  match  closely. 

Once  all  the  one  minute  samples  are  normalized,  they  are  combined  into  a  his¬ 
togram  with  81  bins  spaced  logarithmically  from  0.01  to  100,  or  -40  dB  to  40  dB 
in  power  relative  to  A  —  1.  This  dynamic  range  is  found  to  include  virtually  every 
data  point  and  is  consistent  with  the  range  used  for  previously  published  data. 


3.4  Data  Analysis 

This  section  presents  the  general  results  found  when  fitting  all  three  models  to  all  of 
the  data,  followed  by  specific  results  for  each  measurement  location.  This  is  followed 
by  a  discussion  of  each  of  the  models,  noting  the  range  of  parameters  each  model 
uses  in  fitting  the  data  and  how  these  parameters  vary  with  season,  time  of  day,  and 
station.  Finally,  the  results  are  compared  with  the  voltage  deviation  parameter  Vd. 

Because  of  the  difficulty  of  finding  noise  samples  that  are  completely  free  of  man¬ 
made  interference  and  the  extensive  processing  time  required  to  manipulate  even  the 
uncontaminated  data,  only  a  small  subset  of  the  noise  survey  database  is  analyzed. 
Nonetheless,  the  resulting  data  set  consists  of  35  gigabytes  and  spans  various  times 
of  the  year  and  all  times  of  the  day  at  seven  of  the  stations.  Space  limitations 
preclude  the  presentation  of  the  thousands  of  graphs  and  model  parameter  values 
resulting  from  the  analysis;  instead  the  qualitative  features  of  the  data  are  discussed 
and  some  numerical  results  are  provided  as  proof. 

The  data  clearly  reveal  a  pattern  defining  which  model  works  best  under  which 
conditions,  and  the  findings  are  as  follows:  the  Hall  model  is  found  to  be  very 
accurate  in  modeling  the  amplitude  pdf  of  VLF  radio  noise  under  the  condition 
of  heavy  sferic  activity  in  the  noise,  otherwise  the  a-stable  model  is  optimal.  In 
addition,  there  is  a  fairly  large  transition  region  (as  a  function  of  time,  location, 
frequency,  etc.)  where  both  models  are  equally  accurate.  The  F-L  model  exhibits 
an  error  metric  roughly  10  to  100  times  higher  than  the  other  two  for  most  of  the 
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data  samples  examined,  but  it  is  still  reasonably  accurate. 

As  mentioned  in  Chapter  1,  heavy  sferic  activity  means  that  the  noise  contains 
numerous  sferic  impulses,  as  in  Figure  1.3.  This  condition  is  met  in  general  at 
the  low-latitude  sites  compared  to  the  high-latitude  sites,  as  seen  in  the  differences 
between  the  time-series  data  from  Grafton,  Dunedin,  Stanford,  Kochi  and  L’Aquila 
and  the  data  from  Arrival  Heights,  S0ndrestr0m  and  Thule.  The  higher  latitude 
sites  are  far  from  the  main  sources  of  sferics  (the  tropical  regions)  so  only  larger 
sferics  reach  these  locations.  As  viewed  from  a  high-latitude  receiver,  then,  sferics 
appear  to  be  more  sporadic. 

It  is  stated  in  Section  2.2  that  sferic  activity  at  a  given  location  has  both  a 
seasonal  and  a  diurnal  cycle.  Seasonal  variations  peak  during  the  local  summer,  and 
diurnal  variations  peak  in  the  local  afternoon.  Since  the  optimal  model  depends  on 
sferic  activity,  it  is  thus  related  to  season  and  time  of  day. 

The  sferics  themselves  have  a  frequency  response  (as  seen  at  the  receiver)  that 
peaks  in  the  8-14  kHz  range  and  decreases  with  increasing  frequency;  therefore, 
only  the  stronger  sferics  create  large  impulses  at  the  higher  VLF  frequencies  (as  seen 
in  spectrograms  similar  to  Figures  1.2  and  1.3).  Because  of  this  frequency  response, 
narrowband  time-series  data  appear  to  contain  lower  sferic  activity  as  the  center 
frequency  is  increased. 

The  bandwidth  of  the  receiver  is  known  to  affect  the  impulsiveness  of  the  noise 
because  a  narrow  bandwidth  spreads  the  sferics  in  time,  causing  them  to  overlap 
more  and  appear  less  impulsive  (this  typically  lowers  the  Vd  as  well).  Thus  the 
a-stable  model  exhibits  a  larger  error  than  the  Hall  at  low  bandwidths  since  the 
impulses  are  less  distinctive  against  the  background  noise,  but  this  is  true  only  for 
very  small  bandwidths. 

The  following  section  presents  specific  results  related  to  the  above  discussion, 
but  it  must  be  noted  again  that  the  analysis  is  known  only  to  apply  to  atmospheric 
noise.  Since  auroral  hiss  has  the  characteristics  of  thermal  noise,  it  is  likely  that  the 
a-stable  distribution  with  a  «  2  would  model  it  accurately,  but  this  hypothesis  is 
not  confirmed. 
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3.4.1  Results  by  Station 

The  Arrival  Heights  VLF  data  are  taken  from  the  time  period  May  1995  through 
June  1996;  the  center  frequencies  analyzed  span  the  range  17  -  27.5  kHz  and  the 
noise  bandwidth  is  chosen  to  be  600  Hz  (in  order  to  reject  adjacent  man-made 
signals).  A  fairly  wide  uncontaminated  band  from  25.5  -  27.5  kHz  is  analyzed  as 
well.  It  is  found  that  the  a-stable  model  performs  optimally  over  the  whole  data 
set,  with  an  average  MSLE  of  0.0008,  compared  to  0.006  for  the  Hall  and  0.0153 
for  the  F-L.  This  means  that  on  average  the  a-stable  model  is  roughly  within  7% 
(lO'/o-0008)  of  the  true  pdf,  compared  to  19%  for  the  Hall  and  33%  for  the  F-L.  The 
a-stable  model’s  accuracy  is  especially  good  above  20  kHz,  where  the  respective 
average  errors  are  5%,  22%  and  29%. 

The  Dunedin  VLF  data  cover  the  year  1989;  the  center  frequencies  analyzed 
span  the  range  15  -  27  kHz  and  the  noise  bandwidth  is  400  Hz.  The  band  25  - 
27  kHz  is  analyzed  as  well.  It  is  found  that  the  Hall  model  performs  best  below 
approximately  23  kHz;  above  this  range  the  a-stable  model  is  slightly  better.  The 
respective  percentage  errors  for  the  a-stable,  Hall  and  F-L  pdf’s  are  14%,  10%  and 
33%  for  frequencies  below  23  kHz,  and  8%,  10%  and  31%  for  those  above. 

The  Thule  data  include  June  1986  to  February  1987,  with  center  frequencies 
ranging  from  15  -  43  kHz  (LF  data  were  collected  at  Thule).  The  bandwidth  for  23 
kHz  and  below  is  400  Hz;  at  the  higher  frequencies  tested  it  is  2  kHz.  It  is  found  that 
the  Hall  and  a-stable  models  are  comparable  below  23  kHz  except  during  the  peak  of 
the  seasonal  and  diurnal  cycles,  when  the  Hall  model  exhibits  1/2  to  1/3  the  MSLE 
of  the  a-stable  model.  The  respective  average  percentage  errors  for  the  a-stable, 
Hall  and  F-L  models  below  22  kHz  are  12%,  10%  and  34%;  at  higher  frequencies 
those  errors  are  6%,  15%  and  35%.  Even  at  higher  frequencies  the  seasonal  and 
diurnal  variations  have  an  effect:  at  36  kHz  the  a-stable  error  varies  roughly  from 
5%  to  11%  seasonally  and  from  9%  to  15%  diurnally  (during  the  seasonal  peak). 

The  Spndrestrpm  data  include  September  1993  to  June  1994,  with  frequencies 
from  17  -  26  kHz  and  bandwidths  of  600  Hz.  The  accuracy  of  the  models  follows  the 
same  diurnal,  seasonal  and  frequency  patterns  as  at  the  stations  discussed  above; 
sample  results  for  Spndrestrpm  are  shown  in  Figure  3.2.  The  three  rows  of  this 


Figure  3.2:  Errors  of  the  Hall  (dashed  line)  and  a-stable  (solid  line)  models  as  a 
function  of  frequency  and  season  at  S0ndrestr0m. 
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figure  correspond  to  center  frequencies  of  17.350,  22.550  and  26.500  kHz,  and  the 
four  columns  cover  four  consecutive  seasons  from  September  1993  to  June  1994. 
The  plots  are  of  the  average  percentage  error  as  calculated  using  io VMSLE.  ^g  jja^ 
error  is  depicted  by  the  dashed  line  and  the  a-stable  error  by  the  solid  one. 

Note  that  the  a-stable  model  has  considerably  less  error  than  the  Hall  for  the 
months  of  December  and  March,  when  the  northern  hemisphere  seasonal  variation 
is  at  a  minimum.  However,  in  June  and  September  (near  the  peak  of  storm  season) 
and  at  17  kHz,  the  Hall  model  gives  better  performance.  In  addition,  there  is  a 
diurnal  variation  of  the  a-stable  error  in  September  and  June,  when  it  increases 
with  the  passing  of  North  American  storms  at  approximately  00  UT. 

The  Grafton,  Stanford  and  L’Aquila  data  cover  the  same  general  frequency  and 
time  ranges  as  the  data  presented  above,  with  the  same  results.  Grafton  is  espe¬ 
cially  close  to  heavy  storm  activity  in  North  America,  which  typically  occurs  in  the 
northern  hemisphere  summer,  so  the  accuracy  of  the  a-stable  model  is  noticeably 
worse  during  these  months.  This  is  depicted  in  Figure  3.3,  which  shows  the  average 
percentage  errors  (averaged  over  all  frequencies)  of  the  three  models  as  a  function 
of  time  of  day  for  the  months  of  January  1988  and  July  1988.  The  a-stable  model  is 
clearly  the  best  in  January,  but  in  July  its  performance  degrades  to  a  greater  error 
than  the  Hall  model. 

Bandwidth  Effects  Thule  data  at  center  frequencies  of  36.6  kHz  and  43.4  kHz 
and  with  bandwidths  of  25,  50,  100,  200,  400,  800  and  1600  Hz  were  processed  in 
order  to  test  the  effect  of  increasing  bandwidth.  The  Hall  model  exhibits  a  larger 
MSLE  and  the  a-stable  model  a  smaller  one  as  bandwidth  increases,  but  the  effect 
is  only  seen  during  increased  storm  activity,  and  the  maximum  error  of  either  model 
is  only  approximately  13%. 

An  additional  sample  of  Arrival  Heights  May  1995  data  was  processed  with  a 
center  frequency  of  8  kHz  and  bandwidth  varying  from  25  Hz  to  1600  Hz.  In  this  case 
the  a-stable  model  is  only  7.5  percent  in  error  across  the  whole  range  of  bandwidths, 
while  the  Hall  error  increases  from  9.5  to  21  percent  with  increasing  bandwidth.  The 
F-L  error  varies  between  40  and  60  percent. 
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(a)  January  1988 


(b)  July  1988 


Figure  3.3:  Errors  of  the  Hall,  Field-Lewinstein  and  a-s table  models  as  a  function 
of  frequency  and  season  at  Grafton,  New  Hampshire. 
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Model 

Parameters 

MSdBE 

%  error 

Hall 

ro  =  5.23 

7  =  1.348 

0.0075 

22% 

F-L 

a  =  0.58 
r0  =  0.77 
r  =  0.18 

0.0132 

30% 

a-stable 

a  =  1.844 

7  =  0.293 

0.000038 

1.4% 

Table  3.1:  Parameter  values  for  the  Hall,  Field-Lewinstein  and  a-stable  models  in 
Figure  3.4. 


Model 

Parameters 

MSdBE 

%  error 

Hall 

m  =  3.05 

7  =  0.665 

0.000088 

2.2% 

F-L 

a  =  0.59 
r0  =  0.40 
r  =  0.40 

0.0219 

41% 

a-stable 

o 
"  II 

0.0011 

8% 

Table  3.2:  Parameter  values  for  the  Hall,  Field-Lewinstein  and  a-stable  models  in 
Figure  3.5. 

3.4.2  Parameters  of  the  Models 

This  section  discusses  how  well  each  model  fits  various  portions  of  the  dynamic  range 
of  the  noise  envelope.  This  information  is  not  contained  in  the  MSLE  since  it  is  an 
average  over  the  entire  dynamic  range;  therefore,  this  section  provides  additional 
insight  as  to  why  a  model  may  perform  well  or  poorly.  In  addition,  it  is  stated 
for  each  model  the  range  of  parameters  exhibited  in  fitting  the  data  histograms, 
and  whether  or  not  these  parameters  are  correlated  with  location,  center  frequency, 
bandwidth,  time  of  day  or  season  of  the  year. 

Figure  3.4  shows  the  fit  of  the  three  models  (both  APD  and  pdf)  using  a  typical 
data  sample  for  which  the  a-stable  model  is  optimal;  the  data  and  a-stable  plots 
virtually  coincide.  The  data  are  from  Arrival  Heights  during  the  04  -  08  UT  diurnal 
time  period  in  May,  a  time  and  month  of  relatively  low  sferic  activity.  The  center 


%  exceeding  amplitude  value 
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(a)  fit  of  pdf 


Grafton  APD,  17.5  kHz,  400  Hz  BW,  Mar  1938, 18  UT 


(b)  fit  Of  APD 


Figure  3.5:  Fit  of  the  Hall,  Field-Lewinstein  and  a-stable  pdf  models  to  a  sample 
of  Grafton  data. 
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frequency  is  22.7  kHz  and  the  bandwidth  is  600  Hz.  Parameters  and  errors  for  the 
three  models  are  given  in  Table  3.1. 

The  envelope  pdf  of  Figure  3.4  is  essentially  Rayleigh  except  for  a  heavy  tail  due 
to  occasional  sferics,  starting  at  10  dB  in  the  dynamic  range.  The  Hall  and  F-L 
models  fit  the  curve  accurately  in  the  high  probability  -10  to  10  dB  range,  but  are 
unable  to  accurately  model  the  higher  range.  In  fact,  they  are  often  off  by  orders  of 
magnitude. 

Figure  3.5  shows  the  fit  (both  APD  and  pdf)  for  data  typical  of  when  the  Hall 
model  is  optimal.  The  data  are  from  Grafton,  in  March,  and  contain  heavy  sferic 
activity.  The  center  frequency  is  17.5  kHz  and  the  bandwidth  is  400  Hz.  Parameters 
and  errors  for  the  three  models  are  given  in  Table  3.2. 

The  value  of  m  for  the  Hall  model  and  the  value  of  a  for  the  a-stable  model  are 
significantly  less  in  Table  3.2  than  in  Table  3.1.  The  average  amplitude  is  higher 
relative  to  the  background  noise,  and  so  more  of  the  probability  occurs  at  dB  values 
less  than  zero.  In  addition,  the  data  curve  has  no  convex  bend  at  10  dB  as  in 
Figure  3.4,  making  the  Hall  pdf  an  almost  perfect  fit.  The  upward  curve  of  the 
cc-stable  pdf  near  40  dB  in  Figure  3.5  is  due  to  limits  of  numerical  accuracy  for  the 
algorithm  used  to  determine  the  pdf. 

3.4.3  Error  Across  Dynamic  Range 

The  errors  across  dynamic  range  in  Figures  3.4  and  3.5  are  shown  in  Figure  3.6.  A 
value  of  one  on  the  y-axis  corresponds  to  an  order  of  magnitude  in  error,  so  it  is  seen 
that  the  Hall  and  F-L  models  give  an  estimate  that  is  several  orders  of  magnitude 
too  low  for  the  high  end  of  the  dynamic  range  in  the  Arrival  Heights  data.  The  F-L 
model  gives  too  low  an  estimate  at  low  levels  as  well. 

For  very  small  a,  the  F-L  model  has  functional  form  aa+1,  since  it  is  the  convo¬ 
lution  of  a  with  a“-1  (ignoring  constants).  This  is  clearly  seen  in  Figures  3.4(a)  and 
3.5(a),  where  the  slope  of  the  F-L  curves  at  low  values  is  roughly  0.8  dB/dB.  Since 
the  abscissa  is  in  dB  power,  this  is  equivalent  to  1.6  dB/dB  in  amplitude,  which  is 
almost  exactly  the  value  of  1  +  a.  The  other  curves  (including  the  data)  have  a  slope 
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Arrival  Heights  pdf  error,  22.7  kHz,  600  Hz  BW,  May  1995,  06  UT 


(a)  error  in  Arrival  Heights  pdf  fit 


Grafton  pdf  error,  17.5  kHz,  400  Hz  BW,  Mar  1988, 18  UT 


(b)  error  in  Grafton  pdf  fit 


Figure  3.6:  Errors  of  the  Hall,  Field-Lewinstein  and  a-stable  pdf  models  in  fitting 
the  data  of  Figures  3.4  and  3.5. 
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of  0.5  dB/dB  power,  or  1.0  dB/dB  amplitude,  indicating  a  Rayleigh  pdf  for  small  a. 
At  the  high  end  of  the  dynamic  range,  the  F-L  model  decreases  power-exponentially 
as  a  — y  oo,  which  is  too  rapid  of  a  decay.  The  polynomial  decay  exhibited  by  the 
other  two  models  is  more  accurate. 

3.4.4  Typical  Parameter  Values 

Typical  values  of  the  F-L  parameters  are  as  follows:  a  is  usually  in  the  range  0.5 
to  0.6,  although  values  as  low  as  0.34  and  as  high  as  1.15  are  found.  The  values 
of  r0  and  r  both  range  from  0.1  to  1.0,  but  typically  take  values  of  0.7  and  0.25, 
respectively.  There  is  only  a  mild  dependence  of  the  three  parameters  on  season 
and  time  of  day  at  the  high-latitude  sites,  but  at  lower  latitudes  there  is  a  large 
dependence.  At  Grafton,  for  example,  the  value  of  a  jumps  from  0.5  in  January  to 
0.9  in  June,  and  there  is  an  additional  0.2  diurnal  variation  about  these  values  as 
well.  The  term  r0  correspondingly  drops  from  0.8  to  0.2,  and  r  increases  from  0.2 
to  0.8,  from  January  to  June. 

The  Hall  model  has  values  of  m  that  range  from  2.2  to  11,  with  typical  values 
between  3  and  4.  Values  of  7  range  from  0.6  to  2.4  and  are  typically  1  to  1.5.  There 
is  no  strong  relationship  between  the  Hall  parameters  and  time,  location,  etc.,  but 
the  values  of  m  and  7  do  tend  to  rise  as  the  center  frequency  is  increased  throughout 
the  15  -  27  kHz  range. 

The  a-stable  model  has  values  of  7  that  are  mostly  near  0.3.  The  value  of  a 
is  usually  between  1.6  and  2.0,  but  can  be  as  low  as  1.1.  A  seasonal  and  diurnal 
dependence  is  seen  in  a,  but  it  is  strong  only  at  the  low-latitude  sites. 

3.4.5  ELF  Results 

The  ELF  data  are  severely  contaminated  by  power  line  harmonics  spaced  either 
50  Hz  or  60  Hz  apart,  and  the  need  for  total  rejection  of  these  signals  allows  only 
for  very  narrow  analysis  bandwidths.  (The  need  for  total  rejection  of  power  line 
harmonics  is  demonstrated  by  the  fact  that  even  small  leakage  causes  a  significant 
drop  in  Vj).  Fortunately,  several  data  samples  were  found  for  which  the  240  Hz 
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harmonic  did  not  exist  at  all  for  an  entire  month,  so  the  usable  bandwidth  spanned 
the  range  180  -  300  Hz.  These  data  samples  are  Arrival  Heights  for  September 
1994,  and  Thule  and  S0ndrestr0m  for  April  1990.  They  are  analyzed  using  a  center 
frequency  of  240  Hz  and  bandwidths  ranging  over  the  values  5,  10,  15,  20,  25  and 
30  Hz. 

The  ELF  results  do  not  differ  significantly  from  the  VLF  results  either  qualita¬ 
tively  or  quantitatively.  The  average  percentage  errors  of  the  ar-stable,  Hall  and  F-L 
models  at  Arrival  Heights  are  4%,  12%  and  36%,  respectively;  at  S0ndrestr0m  they 
are  9%,  6%  and  42%  respectively;  at  Thule  they  are  8%,  7%  and  37%.  These  num¬ 
bers  are  averaged  across  bandwidth;  at  the  larger  bandwidths  the  a-stable  model 
always  outperforms  the  Hall.  The  parameter  ranges  of  the  models  are  the  same  as 
described  above  for  VLF. 


3.5  Voltage  Deviation 

The  Vd  statistic  is  reintroduced  at  this  point  because  it  is  a  fundamental  concept 
in  previous  work  on  modeling  atmospheric  noise  [26,  42],  It  is  computed  as  Vd  = 
10  log {E[A?\j E2[A\),  and  as  such  it  is  the  RMS  value  of  the  noise  envelope  divided 
by  the  average  value,  in  dB.  Vd  is  a  measure  of  the  spikiness  of  the  noise  (since 
sharper  impulses  result  in  a  higher  RMS  value  relative  to  the  mean  value)  and  is 
greatly  dependent  on  the  noise  bandwidth;  however,  it  is  not  necessarily  an  indicator 
of  heavy  sferic  activity.  For  the  data  of  Figure  3.4,  the  Vd  is  3.48  dB,  and  for  the 
data  of  Figure  3.5,  it  is  4.91  dB,  but  there  are  many  cases  where  the  o-stable  model 
is  optimal  for  higher  Vd’s  as  well.  (It  should  be  noted  that  the  ability  to  compute 
10  log (E[A2]/ E2[A])  from  the  data  is  based  on  system  dynamic  range  limitations, 
since  Vd  is  undefined  for  many  parameter  values  of  the  Hall  and  a-stable  models). 

Seasonal  and  diurnal  variations  of  Vd  are  seen  in  the  data,  but  they  cannot  be 
systematically  correlated  with  known  storm  distributions.  For  instance,  Vd  does  not 
necessarily  rise  and  fall  with  storm  season:  at  Grafton  (for  one  sample  frequency 
band)  it  reaches  a  low  of  2.5  dB  in  both  January  and  July,  from  approximately  4.5 
dB  in  spring  and  fall.  Other  data  show  that  Grafton’s  Vd  varies  out  of  phase  with  its 
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diurnal  variation,  while  Dunedin’s  Vd  variation  is  in  phase  with  its  diurnal  variation. 
The  typical  range  of  Vd  values  is  1.5  dB  to  6.0  dB,  although  values  as  high  as  10  dB 
can  be  seen  at  high-latitude  sites. 

Vd  is  also  found  to  relate  somewhat  to  the  parameters  of  the  three  pdf  models. 
For  one  set  of  Arrival  Heights  data  samples,  the  Hall  parameters  m  and  7  decrease 
from  10  to  2.8  and  2.0  to  0.4,  respectively,  as  Vd  increases  from  1.5  to  9.0.  The  F-L 
parameters  a,  r0  and  r  change  from  0.7  to  0.4,  1.0  to  0.4,  and  0.1  to  0.2,  respectively, 
as  Vd  increases,  and  the  a-stable  parameters  a  and  7  decrease  from  «  2.0  to  1.2  and 
0.3  to  0.2,  respectively.  The  errors  of  the  three  models  have  fairly  little  relation  to 
Vd,  except  that  the  Hall  model  tends  to  be  more  accurate  than  the  a-stable  model  for 
lower  values  of  m,  which  corresponds  to  higher  values  of  Vd-  This  is  not  universally 
true,  though,  since  as  mentioned  previously  it  is  possible  to  have  a  high  Vd  for  noise 
of  relatively  low  sferic  activity. 

3.6  Conclusions 

This  chapter  presents  many  APD  results  derived  from  a  statistical  analysis  of  low- 
frequency  radio  noise.  Three  models  for  the  noise  envelope  APD  are  described,  and 
the  parameters  and  accuracy  of  each  model  are  determined  as  a  function  of  location, 
time  and  frequency.  The  parameters  and  errors  of  each  model  are  found  to  vary  with 
thunderstorm  activity,  and  the  noise  frequency  and  bandwidth. 

The  main  conclusion  of  this  chapter  is  that  the  Hall  model  is  the  optimal  model 
in  terms  of  accuracy  and  simplicity  for  locations  exposed  to  heavy  sferic  activity,  and 
the  a-stable  model  is  best  for  locations  relatively  distant  from  heavy  sferic  activity. 
A  general  rule  based  on  the  results  would  be  to  use  the  a-stable  pdf  in  high-latitude 
regions  except  at  the  peak  of  the  diurnal  and  seasonal  storm  cycle,  and  to  use  the 
Hall  model  at  low-  and  mid-latitudes  except  at  the  null  of  the  seasonal  and  diurnal 
cycle.  The  crossover  point  is  not  critical;  there  are  a  broad  range  of  conditions  where 
both  models  have  roughly  the  same  accuracy. 


Chapter  4 


Clustering  Poisson  Model 


4.1  Introduction 

This  chapter  presents  a  statistical-physical  radio  noise  model,  denoted  the  “cluster¬ 
ing  Poisson”  model,  and  shows  that  this  model  predicts  several  important  statistical 
features  of  atmospheric  radio  noise  very  accurately.  The  model  is  derived  in  con¬ 
junction  with  the  statistical  results  of  Chapters  2  and  3,  and  so  it  is  supported  by 
a  multitude  of  real  data.  The  term  “clustering  Poisson”  is  used  to  describe  the 
model  because  noise  source  events  occur  as  Poisson  processes  that  are  triggered  by 
another,  independent  Poisson  process,  and  thus  the  noise  impulses  are  seen  to  occur 
in  bursts,  or  clusters. 

A  number  of  statistical-physical  noise  models  (i.e.,  based  on  the  underlying 
physical  process  that  creates  the  noise)  have  been  developed  over  the  past  four 
decades;  each  can  be  roughly  categorized  into  one  of  two  types:  (i)  simple  enough 
to  provide  concise,  closed  form  answers  but  only  somewhat  representative  of  the 
true  physical  situation,  and  (ii)  representative  of  the  true  physical  situation  but 
very  complicated  to  use.  The  second  type  often  uses  specific  storm  distribution  and 
propagation  information  and  calculates  the  noise  characteristics  at  a  given  point 
numerically  ( e.g Warber  and  Field  [46]),  whereas  the  first  type  usually  assumes 
independence  in  space  and  time  of  the  source  distribution  (a  condition  known  not 
to  be  true  in  practice).  The  new  model  partially  bridges  the  gap  between  the  two 
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types  by  replacing  independence  in  time  with  a  clustering  Poisson  assumption. 

The  model  is  verified  by  comparing  it  to  four  statistical  features  of  atmospheric 
radio  noise:  impulse  interarrival  distributions,  correlations  between  impulses,  APDs, 
and  power  spectral  densities.  The  impulse  interarrival  time  distributions  predicted 
by  the  model  are  shown  to  match  those  seen  in  the  data,  thus  verifying  the  accuracy 
of  the  clustering  Poisson  assumption.  Furthermore,  the  noise  amplitude  probability 
density  function  predicted  by  the  clustering  Poisson  model,  the  a-stable  envelope 
pdf,  is  shown  in  Chapter  3  to  accurately  fit  histograms  of  real  data  for  locations 
relatively  distant  from  storm  activity.  Finally,  clustering  of  pulses  is  shown  to  have 
little  effect  on  the  power  spectral  density  of  the  noise,  a  fact  also  seen  in  the  data. 
Given  the  accuracy  to  which  the  predicted  statistical  features  fit  the  actual  data,  it 
is  concluded  that  the  clustering  Poisson  model  is  a  good  candidate  as  a  model  for 
atmospheric  noise. 

Section  4.2  discusses  several  existing  models  for  atmospheric  noise  and  presents 
the  new  clustering  Poisson  model,  then  Section  4.3  shows  that  this  model  predicts 
the  noise  envelope  pdf  to  be  an  o-stable  distribution.  Section  4.4  derives  the  pre¬ 
dicted  impulse  interarrival  distributions,  and  Section  4.5  discusses  power  spectral 
densities.  Finally,  the  results  are  summarized  in  Section  4.6. 


4.2  Review  of  Existing  Models 

Most  impulsive  noise  models  do  not  take  impulse  interarrival  time  dependencies  into 
account,  but  two  statistical-physical  models  that  do  address  impulse  clustering  are 
those  of  Furutsu  and  Ishida  [16]  and  Giordano  and  Haber  [18].  Furutsu  and  Ishida 
address  only  the  clustering  of  pilots  and  leaders  (on  the  order  of  milliseconds)  in 
an  individual  lightning  stroke,  however,  and  the  analysis  does  not  provide  concise 
results.  Giordano  and  Haber  model  impulse  clustering  by  assuming  that  with  each 
independent  impulse  there  is  a  finite  probability  of  a  similar  dependent  impulse  a  time 
Ti  later,  where  r;  is  an  appropriately  specified  random  variable.  The  analysis  is  not 
extended  to  multiple  dependent  impulses  due  to  the  added  mathematical  complexity. 
Both  models  provide  several  key  features  upon  which  the  new  model  is  based,  such 
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as  specification  of  the  source  distribution  and  propagation  characteristics. 

Nikias  and  Shao  [37]  use  an  adapted  form  of  Giordano’s  model  (without  impulse 
clustering),  in  addition  to  a  method  in  Zolotarev  [48]  for  the  model  of  point  sources  of 
influence,  to  prove  that  atmospheric  noise  has  an  a-stable  distribution.  Since  some 
of  the  results  in  the  following  sections  use  this  proof,  the  notation  in  describing  both 
the  non-clustering  model  and  the  new  clustering  model  largely  follows  that  in  [37]. 

4.2.1  A  Non-Clustering  Statistical-Physical  Poisson  Model 

This  section  defines  a  non-clustering  statistical-physical  Poisson  model  that  is  com¬ 
mon  in  principle  to  several  described  in  the  literature  [18].  Begin  by  assuming  that 
the  received  noise  is  the  superposition  of  many  impulses  produced  by  many  sources 
in  a  region  encompassing  the  receiver  location.  Given  this  assumption,  if  the  exact 
source  distribution  were  known,  as  well  as  the  exact  time  that  each  source  emits  an 
impulse  and  the  exact  waveform  of  each  impulse  at  the  receiver  (including  knowl¬ 
edge  of  time  delay),  the  characteristics  of  the  received  noise  waveform  could  be 
completely  determined.  The  problem,  of  course,  is  that  it  could  only  be  determined 
numerically  and  would  not  in  general  yield  simple,  closed  form  results.  In  order 
to  proceed  to  such  results,  exact  knowledge  of  the  source  characteristics  must  be 
replaced  with  simplified  statistical  approximations  and  expectations. 

The  region  Q  encompassing  the  receiver  location  is  defined  on  Rn,  where  n  € 
{1,2,3}  is  the  dimension  of  the  space  ( i.e a  line  (R1),  a  plane  (R2),  or  a  volume 
(R3)).  The  receiver  is  at  the  origin  and  each  source  i  is  at  a  position  x,  a  distance 
|x;|  from  the  receiver.  It  is  assumed  that  all  sources  have  similar  enough  waveform 
generation  mechanisms  that  their  emitted  waveforms  may  be  modeled  as  a,\D(f ;  #;), 
where  the  random  amplitudes  a;  are  independent  and  identically  distributed  (i.i.d.) 
with  pdf  fa(a)  (denoted  i.i.d.  ~  fa(a))  and  the  random  parameters  0,-  are  i.i.d.  ~ 
/©(0).  The  0i  represent  arbitrary  mappings  from  a  probability  space  to  an  ensemble 
of  possible  waveforms;  they  are  chosen  to  be  scalars  for  simplicity. 

The  effect  of  the  transmission  medium  is  modeled  as  the  combination  of  a  power 
law  attenuation  factor  u  and  a  linear,  time-invariant  (LTI)  filtering  factor  h(0)  such 


CHAPTER  4.  CLUSTERING  POISSON  MODEL 


49 


that  the  impulse  a;Z)(t; 0t)  appears  as  9i)  *  h(9{)  at  the  receiver,  where  ct 

is  a  positive  constant  and  u>  0.  Since  0;  is  arbitrary,  the  convolution  D(t;  9i)*h(9i ) 
may  be  defined  as  E(t;  8t)  without  loss  of  generality.  Further  defining  the  quantity 

0i)  =  t \E(t]  9i )  (4.1) 

lX»l 

and  letting  the  random  variable  N  be  the  number  of  contributing  impulses  at  the 
time  of  observation,  the  received  waveform  Y  is  then 

N 

y  =  X) aiU  (*«•;  &i)-  (4.2) 

*=1 

The  time  >  0  is  defined  as  the  difference  between  the  observation  time  t  =  0  and 
the  source  emission  time,  so  it  is  larger  for  older  impulses  ( i.e t  is  a  measure  of 
negative  time). 

It  now  remains  to  specify  N.  In  order  to  attain  the  most  tractable  results  while 
maintaining  a  reasonably  accurate  physical  representation,  N  is  assumed  to  be  the 
number  of  events  generated  by  a  Poisson  process  in  space  and  time  with  source 
density  function  p(x,t).  The  value  p(x,t)dxdt  is  then  the  probability  that  a  noise 
impulse  will  be  emitted  from  the  (infinitely  small)  region  defined  by  the  line,  square 
or  cube  (for  n  =  1,2,3)  with  far  corners  x  and  x  +  dx,  and  during  the  (infinitely 
small)  interval  [t,t  +  dt]  prior  to  observation.  The  term  p(x,  t)  may  take  a  general 
form  over  x  €  Cl  and  t  6  [0,  oo),  but  for  simplicity  it  is  approximated  as 

(4.3) 

where  po,fi  >  0  are  constants.  Thus  the  sources  are  defined  to  be  independent  both 
in  time  and  direction  from  the  receiver.  The  exponent  p  defines  a  variation  in  source 
density  with  distance,  providing  an  added  degree  of  freedom  with  minimal  added 
complexity. 

The  definition  of  p(x,  t )  in  Eq.  (4.3)  is  the  key  approximation  for  making  further 
analyses  reasonably  simple  and  tractable,  but  it  is  not  the  true  physical  source 
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distribution  of  atmospheric  radio  noise.  Various  storm  centers  are  active  at  various 
times  of  day  and  year  [4,  5,  19],  so  for  a  given  receiver  location  there  are  certain 
directions  from  which  heavier  sferic  activity  arrives  than  others.  This  directional 
dependence  contradicts  the  source  distribution  of  Eq.  (4.3),  but  nonetheless  Eq.  (4.3) 
is  often  reasonably  accurate. 

The  temporal  independence  approximation,  however,  is  not  reasonably  accurate. 
Bursts  of  sferics  on  the  order  of  one  second  in  length  are  clearly  seen  in  spectrograms 
of  ELF/VLF/LF  time-series  data  at  all  locations  and  times  (see  Figure  1.2),  and  they 
have  considerable  effect  on  the  time-series  data.  The  model  presented  below  adds 
impulse  clustering  properties  to  the  model  just  described  in  a  way  that  accurately 
represents  this  clustering  of  sferics. 


4.2.2  Definition  of  the  Clustering  Poisson  Model 

Instead  of  specifying  the  received  waveform  Y  to  be  the  superposition  of  impulses 
distributed  throughout  the  region,  specify  it  to  be  the  superposition  of  clusters 
distributed  throughout  the  region,  where  a  cluster  is  defined  as  a  burst  of  impulses 
from  a  given  location  over  a  short  time  period.  (Clusters  will  be  specified  in  detail 
in  the  next  two  sections).  The  resulting  waveform  at  the  receiver  is  then 

Nc  Nc  Nk 

Y  =  '£xt  =  ’£’£**, (4.4) 

k—  1  A;=l  2=0 

where  Nc  is  the  number  of  clusters,  Nk  is  the  number  of  impulses  in  cluster  k,  and 
the  indices  k,i  on  a,  r  and  9  refer  to  the  ith  impulse  of  the  Arth  cluster  (i  —  0 
refers  to  the  impulse  that  starts  the  cluster).  All  the  impulses  in  cluster  k  are 
modeled  as  occurring  at  the  same  location  x^.  The  cik.i  and  6k,i  are  all  independent, 
an  assumption  verified  by  the  lack  of  correlations  between  impulses  discussed  in 
Section  2.5.  Every  cluster  is  assumed  to  be  independent  of  every  other  in  both 
space  and  time.  The  X^s  are  the  waveforms  of  the  individual  clusters. 
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4.2.3  Specification  of  Clusters 

The  number  of  clusters  Nc  is  analogous  to  N  in  Eq.  (4.2)  and  the  clusters  have 
the  source  density  of  Eq.  (4.3),  but  the  statistical  properties  of  the  clusters  them¬ 
selves  still  need  to  be  specified.  The  specification  chosen  must  accurately  model  the 
physical  mechanism  of  lightning  and  should  also  lead  to  tractable  results. 

In  order  to  develop  a  physically  accurate  model,  several  features  of  lightning  must 
be  considered.  From  Uman’s  classic  text  Lightning  [44],  it  is  known  that  a  complete 
cloud  to  ground  discharge,  called  a  flash ,  is  made  up  of  one  or  more  intermittent, 
partial  discharges,  called  strokes.  Each  stroke  involves  a  complex  process  whereby 
either  a  stepped  leader  or  a  dart  leader  initiates  a  strong  return  stroke,  transferring 
large  amounts  of  charge  between  the  earth  and  its  atmosphere.  A  significant  per¬ 
centage  of  flashes  contain  many  strokes;  one  data  sample  in  [44]  indicates  that  40 
percent  of  all  flashes  contain  at  least  five  strokes  and  some  contain  as  many  as  19. 
In  addition,  50  percent  of  those  flashes  with  five  or  more  strokes  have  a  duration  of 
more  than  400  milliseconds.  Data  histograms  in  [44]  and  [46]  showing  typical  flash 
per  stroke  distributions  indicate  that  the  number  of  flashes  per  stroke  can  be  rea¬ 
sonably  modeled  as  a  geometric  random  variable;  thus  the  number  of  impulses  per 
cluster  in  the  clustering  Poisson  model  is  defined  to  be  a  geometric  random  variable 
as  well. 

Now  the  timing  of  impulses  within  a  cluster  must  be  specified.  Every  cluster 
has  at  least  one  impulse,  which  is  the  start  of  the  cluster.  The  interarrival  times 
Ti  between  any  additional  impulses  in  a  given  cluster  are  specified  as  i.i.d.  random 
variables  with  pdf 

h{t)  =  A1e-Al<,  t  >  0,  (4.5) 

i.e.,  an  exponential  distribution  (denoted  EXP(Ai)).  Each  cluster  then  is  essen¬ 
tially  a  variable  length  Poisson  process  with  rate  Ai  impulses  per  second,  and  as 
stated  above,  all  the  clusters  are  independent  of  the  main  process  and  of  each  other. 
Similar  cluster  specifications  have  found  use  in  modeling  computer  failures  [29], 
earthquakes  [45],  and  neural  impulses  [21]. 
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4.2.4  Length  of  Clusters 

This  section  specifies  the  probability  mass  function  (pmf)  of  the  NC s  and  determines 
the  pdf  of  the  cluster  lengths  Tzfc-  The  NC s  are  modeled  as  i.i.d.  geometric  random 
variables  with  pmf 


P(N  =  n)  = 


Az 

Aj  +  A  £, 


Ai  +  A  l) 


n  =  0, 1,2, 3 , 


(4.6) 


where  the  index  k  is  omitted  since  this  pdf  applies  for  all  k.  The  parameter  Az  is 
now  shown  to  be  the  reciprocal  of  the  expected  cluster  length. 

If  N  =  0,  there  are  no  impulses  after  the  first  impulse  and  the  cluster  length  is 
zero.  For  a  given  N  >  0,  the  length  Tl  of  the  cluster  is  the  sum  of  N  independent 
EXP(Ai)  random  variables,  and  therefore  it  is  known  to  have  the  Erlang  density 


hL\N{t,n) 


Ai(A 

(n  -  1)! 


Thus  the  pdf  frL(t)  for  t 

hL(t)  = 


>  0  can  be  determined  by  conditioning  on  N: 


yv  Az,  (  At  V  Ai(A1t)n-1e~Al* 
n=i  Ai  +  A l  \Ai  +  Xl)  (n  —  1)! 


(4.7) 


(4.8) 

(4.9) 
(4.10) 


where  (4.10)  is  obtained  using  the  exponential  summation  formula.  Note  that  the 
rest  of  the  probability  on  Tl  is  at  the  point  P(Tl  =  0)  =  A l/{Xi  +  XL). 

Since  Ai/(Ai  +  Az)  is  the  probability  that  N  >  0,  it  is  apparent  from  Eq.  (4.10) 
that  hL\N>o  is  EXP(AiAz/(Ai  +  Az))  distributed.  In  addition,  since  the  expected 
value  of  an  EXP(A)  random  variable  is  1/A,  it  follows  from  (4.10)  that 


E[Tl]  = 


(4.11) 
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An  equivalent  cluster  specification  results  by  defining  clusters  to  be  Poisson 
processes  of  random  length  T^,  where  T^  is  EXP(A.l).  In  other  words,  it  is  equivalent 
to  define  (i)  the  A^’s  as  geometric  random  variables,  so  that  the  cluster  length  Tl 
given  Nk  is  Erlang  distributed  as  in  (4.7),  or  (ii)  the  length  T^  of  a  rate  Ai  Poisson 
process  such  that  Nk\T^  is  a  Poisson^TjJ  random  variable.  With  definition  (ii), 
the  pmf  of  Nk  is 


P(Nk  =  n ) 


Jo  n\ 

-7  r^tre-^dt 

n\  Jo 


A  L  (  Ai  A 
Ai  +  Al  \Ai  +  Ax,y 


n  =  0, 1,2,3..., 


(4.12) 

(4.13) 

(4.14) 


where  the  integral  in  (4.13)  is  solved  using  item  (860.07)  in  Dwight  [9].  Since 
(4.14)  and  (4.6)  are  identical  and  the  interaxrival  specifications  are  the  same  as  well, 
definitions  (i)  and  (ii)  are  equivalent.  The  variable  Tt,  however,  is  not  the  length  of 
the  cluster;  it  is  simply  the  length  of  an  underlying  Poisson  process.  The  variable  Tl 
(which  is  less  than  or  equal  to  X^)  is  the  true  length  of  the  cluster,  for  it  indicates 
the  location  of  the  last  impulse.  Definition  (ii)  is  introduced  because  it  simplifies 
the  proof  in  the  next  section. 


4.3  Probability  Density  Function  of  Clustered  Pois¬ 
son  Noise 

Now  that  the  clustering  Poisson  model  is  specified,  it  remains  to  demonstrate  its 
validity  in  the  modeling  of  real  data.  This  section  presents  a  proof  that  the  first 
order  pdf  of  the  received  noise  Y  defined  by  the  clustering  model  of  Eq.  (4.4)  is  an 
a-stable  distribution. 

Let  the  receiver’s  response  due  to  a  single  cluster  Xk  starting  at  time  tk  before 
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zero  (zero  is  the  observation  time)  and  at  location  Xk  be 

Nk 

Xk  =  dkfiU (tk'i  Xk,  @k,o)  +  ^  ^  ak,iU Xfc,  @k.i)i  (4.15) 

i=  1 

where  the  first  term  of  (4.15)  is  the  impulse  at  4,  the  cluster  starting  time.  Since 
the  clusters  are  i.i.d.,  the  index  k  is  redundant  when  specifying  any  one  cluster;  thus 
the  notation 

N 

X  =  aoU(t]x,0o) +  '^r,aiU(Ti‘,x,0i)  (4.16) 

«=1 

is  used  from  this  point  on,  where  i  =  0, . . . ,  N  remains  the  index  for  impulses  within 
a  cluster.  (In  Eqs.  (4.15)  and  (4.16),  the  i  =  0  term  is  separated  out  from  the 
summation). 

Now  use  the  second  cluster  specification  of  Section  4.2.4,  noting  that  the  length 
1 1  is  an  EXP(Al)  random  variable.  The  impulse  rate  within  clusters  over  this  length 
is  Ai  impulses  per  second,  so  it  follows  that  for  a  given  ti,  N  ~  Poisson(A1<£i)  and 
the  tC s  (for  i  >  0  and  without  respect  to  ordering)  are  i.i.d.  uniformly  distributed 
over  [t  —  <£,t].  (Once  again,  time  is  in  the  negative  direction).  Defining  = 
a,iU{Ti]  x,  as  the  contribution  to  the  received  waveform  of  any  one  of  the  impulses 
other  than  the  first  in  a  cluster  X ,  it  also  follows  that  the  characteristic  function  of 
the  conditional  pdf  of  given  x,  t  and  denoted  $,-(o;)|x,t,  is 

$i(u)\x,t,tL  =  E[eiwa'u{-Ti'*A)]  (4.17) 

=  J^fa(a)da  fQfe(0)d6  J* *  j- dr  ejuaU^e\  (4.18) 

where  the  integrand  in  (4.18)  is  brought  to  the  far  right  for  reading  clarity.  (This 
format  is  used  in  several  of  the  following  equations  as  well).  The  impulsive  waveform 
U  (r;  x,  6)  is  specified  without  loss  of  generality  to  be  zero  for  r  <  0  in  order  that 
the  system  remains  causal. 


The  characteristic  function  of  X  given  x  and  f ,  including  the  i  =  0  (first)  impulse, 
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is  then 


(4.19) 

=  fa(a)da  j&  fe($)d6  e^aV^ 

roo  ,  00  ( \.r)nP~XlT 

■  /  S-!^L - (*.'(u,)|x,t,T)”  dr 

Jo  £t0  n ! 

(4.20) 

—  j  j°°  ^ie-AireAir(($i(a;)|x,f,r)-X) 

(4.21) 

—  f  f  ejuaU(t;x,6) 

.Ja  Je 

■  f°°  Me~XLTel  Ma)daf°  Me)defir 

Jo 

)  dr. 

(4.22) 


Note  that  the  —1  in  the  last  exponent  of  (4.21)  is  brought  inside  the  integrals  over 
a,  6  and  t'  in  (4.22)  (the  joint  pdf  will  still  integrate  to  the  same  result),  and  the 
notation  of  the  bracketed  term  in  (4.20)  is  simply  abbreviated  in  (4.21)  and  (4.22). 

In  order  to  get  the  characteristic  function  of  any  cluster  that  is  randomly  dis¬ 
tributed  over  a  region  x  €  0(Ri,  R2)  (i-e.,  the  region  for  which  R\  <  x  <  R2)  and 
a  starting  time  t  €  [0,  T],  expectations  must  be  taken  over  x  and  t: 


roo 

/  ALe 
Jo 


-A  Lr 


dr>  chtdt,  (4.23) 


where  again  the  integral  expressions  on  a  and  0  have  been  abbreviated. 

Now  define  the  superposition  of  all  clusters  within  the  region  fi(i?i,i?2)  and 
starting  in  the  time  interval  [0,  T]  as 


Nc  Nc  Nk 

Yt,RuR2  )  J  -Affc  —  y  \  y  ]  Gk,iU{'Tk,i]  Xfc,  0k,i) 

k=l  k~l  t=0 


(4.24) 


as  in  Eq.  (4.4).  The  number  of  clusters  Nc  is  a  Poisson  random  variable  with 
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parameter 

=  Lto  D  ,  f  P^)  dtdx  (4.25) 

and  the  (x*,ifc)  terms  (the  locations  and  starting  times  of  the  clusters)  are  i.i.d. 
with  distribution 

fT,R„R,  =  -fhl *1,  X  €  n(R„R,),  t  €  [0,  T).  (4.26) 

Since  cluster  formation  is  a  Poisson  process,  the  characteristic  function  of  Yt,rur2 
can  be  calculated  in  the  same  manner  as  the  second  term  on  the  right  in  (4.19)  to 
get 

$K(w)  =  (4.27) 

Inserting  (4.26)  into  (4.23)  and  the  result  into  (4.27),  then  taking  logarithms,  results 
in 


IM)  =  l  i/M  U/P"]  />e 


=-Atr 


dr 


— 1>  dxdt,  (4.28) 


where  the  last  -1  is  brought  inside  the  integrals  over  x  and  t  without  changing  the 
result.  Inserting  (4.1)  and  (4.3),  noting  that  fa(a )  is  a  symmetric  distribution,  and 
using  the  identity  cos  2x  =  1  —  2sin2(x)  results  in 


log*y(u,)  =  //Xi^|[2/”/eCOs(i,<.o1|x|---£(<;«))][/o00AIe-^ 


(4.29) 


At  this  point  the  last  exponential  in  (4.29)  must  be  converted  to  a  Taylor  series, 
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resulting  in 


log$y(a>)  =  j~  \Le 

(~4Xl)1  dr 


b=o 


—1  >  dxdt. 


(4.30) 


In  order  to  proceed  further,  solve  for  only  the  j  =  0  term  and  group  the  last  —  1 
with  this  term  as  well.  Since  the  integral  over  r  is  one  in  this  case,  the  j  =  0  term 
of  is 


$y,j=o(w) 


=  Jl  In  wt[“4Co  J©  sm2(^aCl|x|  v E(P,8))\  dx  dt . 


(4.31) 


Now  integrate  over  x  first,  and  consider  the  substitution  r  =  |x|.  If  n  =  1,  the 
integral  is  over  dr;  if  n  =  2,  switching  to  polar  coordinates  results  in  an  integral 
over  2nr  dr ;  if  n  =  3,  spherical  coordinates  give  an  integration  over  4irr2  dr.  For 
general  dimension  n,  the  integral  is  over  cnrn~1  dr.  Using  this  conversion  and  an 
additional  substitution  r'  =  r~u,  dr'  =  —ur~u~ldr  (and  then  dropping  the  primes) 
results  in 


_  [/0T  CofefX  s^a^rWmr-^-^drdedadt] 

vY,j= o  —  e  L  J. 


(4.32) 


Now  set  R\  — >  0,  i?2  — *■  oo  and  T  — >  oo,  define  a  =  where  0  <  a  <  2,  and  use 
the  definite  integral  (3.823)  from  Gradshteyn  and  Rhyzik  [20], 


T(fx)  cos  1y 

2^+1  Q/4  ’ 


a  >  0,  —2</j,<0, 


(4.33) 


d=o 


_  g— 7o  M 


to  get  the  result 


(4.34) 
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where 


7o  = 


2c“cnp0r(l  —  a) cos 


7TQ? 

2 


z/a 

roo  /»  /»oo 

/  /  /  aa|E(*;0)|V.(«)/*(0)^<k.  (4.35) 

JO  J©  JO 


(Note  that  the  relation  T(1  —  a)  =  — aT(— a)  is  used  in  (4.35)). 

In  order  to  solve  for  the  j  =  1  term  of  Eq.  (4.30),  rewrite  it  using  this  identity: 


sin2(o)  cos(2 b)  =  sin2(a  +  b)  +  sin2(a  —  b)  —  sin2(5) 

z  z 


(4.36) 


to  get 


log  $y,j=i  (w)  =  r/^-rxLe-^f.8Xlr  f  r  f  f 

JO  J$l  |x|^  Jo  I  Ai  =0  -/©I  i 02=0  i@2  </ T;=f-T 

[^sin2  Qtuci|x|_I/(aiE(t;0i)  +  a2E{r'\  62))j 
+^sin2  {^ci\x\~'/(alE{t\61)  -  a2£(r';  02))^) 


-  sin2  Qwcx|xj  l/alE(t‘,0i)^ 


dr' ...  drdx dt. 


(4.37) 


The  three  sin2()  terms  in  (4.37)  may  be  integrated  separately,  each  in  a  manner 
analogous  to  Eqs.  (4.31)  -  (4.35),  to  determine  the  result 


*rj=  i  =  e-7,|uj|a, 


(4.38) 


where 


7i  = 


4c“c„/)oAir(l  -  a)  cos 


7 ra 
2 


i/a 


r  r^-XLT  r  t  r  f  i* 

JO  Jo  J a\~ 0  J©j  Ja 2=0  */© 2  Jt—T 


CL l  =0  J ©1  JcZ2=0  J ©2 

^ki£(<;  *i)  +  a^r';  *2)|“  +  0i)  -  a2£(r';  02)|“ 


dr'  d6 2  da2  da*  dr  dt. 


(4.39) 
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Moving  on  to  j  =  2,  the  bracketed  term 

L  L  L  dr'  (si“2(5“acilxr‘'£(T'; *)))  <4-40) 

from  Eq.  (4.30)  is  now  squared,  resulting  in 

roc  p  pt  poo  p  pt 

Jc i2=0  j@ 2  Jt-T  2  Ja3=  0  J&3  Jt-T  3 

(sin2(iajo2ci|x|~l/£’(r2;02)))  (sin2(iwa3Ci|x|-1/£(r3;  03)))  .  (4.41) 

Using  the  formula 

111  1 

sin2  a  sin2  b  =  -  sin2  a  +  -  sin2  b  —  —  sin2(a  +  b)  —  -  sin2(a  —  6),  (4.42) 

results  in  an  expression  of  the  form  of  (4.37)  but  now  there  are  nine  sin2()  terms 
instead  of  three.  When  each  is  integrated  as  in  (4.32)  -  (4.35),  the  result  is 


®Y,J=2  =  e-^° 


(4.43) 


where 

72  _  8c°cnp0Af r(l  -  a)  cos  r°°  ^&_XlT 

i/a  Jo  Jo 

roo  p  poo  p  pi  poo  p  pt 

/  fadai  /  fedOi  /  fada2  /  fedd2  /  dr2  /  fada3  /  fed03  /  dr. ' 

Jo  Je-i  Jo  Je 2  Jt-r  Jo  Je3  Jt~r 

“I aiEftOi)  +  a2E(r'-e2) |“  -  \\ a1E(t]dl)  -  a2E(r^62)\a 
-\\ <hE(t;6x)  +  a3E(r'-,e3)\a  -  i| .«,£(*;  0j)  -  a3E(r^d3)\a 
+-  \aiE(t;  0i)  +  a2E(r2; 02)  +  a3E{r^,  #3)!“ 

+  g \a\E(t\  $i)  +  a2E(r2; 62)  —  a3E{r 3;  03)|“ 
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+  g  &i )  —  a,2E(T2, 02)  +  a2,E{r 3;  ^3)!°' 

+  0  0i)  —  &2-£'(r2)  ^2)  —  a3-E’('r3j  ^3)!°} 


Jr  dt. 


(4.44) 


For  the  j  =  3, 4, 5 . . .  terms,  each  increment  of  jr  adds  three  additional  integrations 
and  multiplies  the  number  of  sin2()  terms  (as  in  (4.37))  by  three;  however,  the 
relation  (4.42)  may  be  applied  repeatedly  until  3J  individual  sin2()  terms  remain, 
with  arguments  consisting  of  the  sums  and  differences  of  the  aj£(r,-;  6{)  terms.  It 
follows  that  V  integrations  of  the  form  (4.32)  -  (4.35)  may  then  be  performed  until 
the  result 

=  e~vMa  (4.45) 

is  obtained,  but  the  7 j  become  increasingly  complicated  to  express.  Convergence  is 
assured  since  the  exponential  Taylor  series  converges  for  all  real  arguments  and  by 
inspection  of  (4.29),  from  which  it  is  apparent  due  to  bounds  on  the  sine  and  cosine 
functions  that  the  term  in  braces  is  between  -1  and  0  for  all  x,  t. 

The  final  desired  result  (^)  is  found  by  multiplying 

*  y  (w)  =  n  *rj i=  1 M° .  (4.46) 

j=0 

and  thus  it  is  proved  that  the  characteristic  function  of  the  received  noise  pdf  is 
that  of  an  ct-stable  distribution.  By  the  uniqueness  of  characteristic  functions,  this 
proves  that  the  marginal  pdf  of  the  noise  is  an  a-stable  distribution. 


4.3.1  Probability  Density  Function  of  the  Noise  Envelope 

The  previous  sections  derive  the  pdf  of  the  received  noise,  but  it  remains  to  determine 
the  statistics  of  the  narrowband  envelope.  Since  the  E(t,0)’s  of  Eq.  (4.1)  may  be 
assumed  to  include  an  arbitrary  bandpass  filter  response  without  loss  of  generality, 
it  remains  only  to  prove  that  the  equivalent  lowpass  signals  (as  in  Eq.  (2.1))  are 
jointly  a-stable. 
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Let  E(t;  6)  be  narrowband,  so  that  it  can  be  represented  as 

E(t]0)  =  Ei(P,9)cos(2T:fot)  —  Eq (t;  0)  sm(2n  fat)  (4-47) 

=  A(t ;  0)  cos  <f> cos(2nfot)  —  A(t;  0 )  sin  4>  sin(27r/of ),  (4.48) 

where  A  and  cj)  are  the  complex  amplitude  and  phase  as  in  Section  2.3.  Equation  (4.4) 
may  then  be  written 

Y  =  Yj  cos(27r /0f)  —  Yq  sin(27r  f0t)  (4.49) 

Nc 

=  X  {Xn  cos(27r/0t)  -  XQk  sin(27r/0t)}  (4:50) 

fc=i 
Nc  Nk 

=  EEtS1  {A 9k,i )  cos  <f>k,i  cos(27T  f0t) 
k= 1 i= o  lxl 

-  A(tk/,  0k, i)  sin  <f>k,i  sin(27r/0t)}  .  (4.51) 

Now  assume  the  </>’s  are  i.i.d.  uniform  [0,  27t),  as  in  Section  2.3.  The  joint  char¬ 
acteristic  functions  Qx^Xq^)  and  $yJtyQ(w)  (where  again  u  =  (lji,uj2))  may  be 
determined  using  an  analysis  similar  to  that  of  the  previous  section,  the  equivalent 
of  Eq.  (4.29)  being 

iog*w,)=/;/^ 

[2/  o/e^i  cos  (aci\x.\~v A(t\ 9)(uji  cos<j>  +  u2  sin<^>))  J  A 

g  [_4Al  /o=0  J©  H-T  dr'&fo*  d<t>  sin2(|  aci  cos  {f>+u>2  sin  ^))j 

(4.52) 

Using  a  Taylor  series  as  before  and  noting  that  cos  <p+u2  sin  (f>  may  be  written 

|^|  cos  {(j)  —  4>u)  (where  ^  =  arctan(^)),  the  equivalent  of  Eq.  (4.30)  is 


CHAPTER  4.  CLUSTERING  POISSON  MODEL 


62 


Compaxing  (4.53)  to  (4.30),  it  is  apparent  that  4>y7jyg  may  be  solved  in  the 
manner  of  Eqs.  (4.30)  -  (4.46)  to  determine 

Sy„yM  =  ft  ’’)y\  (4.54) 

j= 1 

where  for  example 


2c°1cnPor(l  —  a)  cos 


poo  p  poo 

/  //  aa\A(i~,e)\afa(a)M6)dtdeda.  (4.55) 

Jo  Je  Jo 


The  characteristic  function  Eq.  (4.54)  defines  an  isotropic  a-stable  random  vec¬ 
tor.  It  is  shown  in  [37]  that  such  a  random  vector  has  a  complex  envelope  A  = 
\jXj  +  Xq  with  the  pdf  of  Eq.  (3.5);  therefore,  it  is  finally  proved  that  the  cluster¬ 
ing  Poisson  model  predicts  an  a-stable  envelope  distribution  for  A. 


4.3.2  Validity  of  Parameters 

The  parameter  v  in  Eq.  (4.1)  is  the  rate  of  attenuation  of  the  noise  signals  with 
distance;  typical  values  used  in  previous  models  range  from  0.5  to  1.0  [18].  Given 
that  n  =  2  is  the  most  realistic  choice  from  a  physical  standpoint  (or  perhaps 
slightly  larger  than  two  for  an  empirical  approximation),  and  fi  is  fairly  small  since 
the  source  density  does  not  decay  significantly  with  distance,  it  follows  that  values 
of  a  in  the  range  1.5  -  2.0  are  quite  reasonable.  Note  that  for  n  >  u  +  p,  the  noise 


CHAPTER  4.  CLUSTERING  POISSON  MODEL 


63 


is  infinite  when  R2  -»  oo,  but  this  is  compensated  for  by  fixing  R2  at  some  large 
value  f?max  with  the  assumption  that  it  does  not  appreciably  affect  the  integrals  of 
the  form  of  (4.32)  -  (4.35). 


4.3.3  Arbitrary  Receiver  Gain  Pattern 


The  preceding  analyses  assume  that  the  receiver’s  gain  pattern  is  omnidirectional, 
i.e.,  the  same  in  all  directions.  It  is  now  shown  that  if  the  receiver  has  an  arbi¬ 
trary  gain  pattern,  the  received  noise  waveform  is  still  a-stable  distributed.  This 
result  follows  because  the  key  integration  that  results  in  the  a-stable  distribution, 
Eqs.  (4.32)  -  (4.35),  is  over  distance  r  from  the  receiver  and  not  over  angle  of  arrival. 

An  arbitrary  gain  pattern  may  be  incorporated  without  loss  of  generality  as  a 
functional  dependence  of  fa(c t)  on  direction  angle  (  (n  =  2  is  assumed,  again  without 
loss  of  generality),  so  that  fa(a)  is  now  /a(c(a,  ().  Using  this  modification,  Eq.  (4.32) 
becomes,  for  example, 


$Y,j= o  =  e 


4cnPQ 


f  fT  r2’ 

IJ t=zO  27 r  Jo 


d<  Co  A|c(».C)  fe  CRl  s*n2(  juiaci  rE(t;6))r  a~ldr 


(4.56) 


which  still  results  in  an  expression  of  the  form 


$yj=o  =  e 


-to  kl“ 

? 


(4.57) 


but  now  70  is 

2c“c„por(l  —  a)  cos  Zf- 
7o  =  - 

VOL 

1  r2TT  roo  r  roo 

'S7o  Jo  J&  Jo  ^W^TM^OfsWdtdedadC  (4.58) 

All  the  other  terms  will  have  the  same  modification  to  their  7 j  terms,  and  the 

final  result  will  still  have  a  characteristic  function  of  the  form  of  (4.46)  and  (4.54). 
Thus  an  arbitrary  gain  pattern  may  be  added  to  the  analysis  without  changing  the 
result  that  the  noise  is  a-stable  distributed. 
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4.4  Impulse  Interarrival  Distribution  of  Clustered 
Poisson  Noise 

This  section  shows  that  the  clustering  Poisson  model  accurately  predicts  the  im¬ 
pulse  interarrival  time  distributions  of  low-frequency  noise.  Predicted  interarrival 
distributions  are  determined  from  the  model,  and  these  are  shown  to  match  the 
distributions  derived  from  the  data. 

As  defined  in  Section  4.2.3,  the  clustering  model  proposes  that  the  number  of 
follow-on  impulses  (after  the  first)  within  a  given  cluster  is  a  geometric  random  vari¬ 
able,  and  that  the  interarrival  times  between  adjacent  impulses  within  each  cluster 
are  independent  and  EXP(Ai)  distributed.  The  total  process  is  the  superposition  of 
all  the  impulses  in  all  the  clusters,  and  the  clusters  themselves  are  occurrences  of  a 
Poisson  random  process  with  parameter  A2  clusters  per  second. 

Since  the  occurrence  of  an  impulse  is  defined  as  the  crossing  of  a  threshold  level 
as  in  Figure  2.7,  the  individual  impulses  can  be  represented  for  simplicity  as  Dirac 
delta  functions.  The  total  process  can  then  be  expressed  as 

00  Nk 

y(i)=  X  X  *(*  -  tk  -  n,i),  (4.59) 

k=— 00  t=0 

where  Nk  is  the  number  of  impulses  in  cluster  k ,  tk  is  the  start  time  of  cluster  Ar, 
and  Tkj  is  the  time  from  the  start  of  the  cluster  k  to  the  ?th  impulse  of  cluster  k. 
(Time  is  now  in  the  forward  direction).  Note  that  Tk, 0  =  0  for  all  k ,  and  also  that 
Nk  —  0  for  clusters  with  no  follow-on  impulses. 

If  a  stopwatch  is  started  at  an  arbitrary  point  in  time  and  stopped  at  the  next 
occurrence  of  an  impulse,  the  time  shown  on  the  stopwatch  is  the  waiting  time  for 
the  next  impulse,  or  the  forward  recurrence  time.  The  forward  recurrence  time  is  a 
random  variable  Tw  with  marginal  pdf  frw{t). 

If  the  stopwatch  is  instead  started  at  the  occurrence  of  an  arbitrary  impulse  and 
stopped  at  the  next  impulse,  the  time  on  the  stopwatch  is  the  impulse  interarrival 
time.  This  random  variable  is  denoted  7/,  with  marginal  pdf  /y7(t).  The  following 
sections  derive  both  fow{t)  and  /r7(t)  for  the  clustering  Poisson  model,  and  then 
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compare  the  derived  fTw{t)  to  the  interarrival  distributions  seen  in  the  data. 


4.4.1  Waiting  Time  pdf 

The  pdf  frw{L)  of  the  forward  recurrence  time  of  the  process  y(t)  is  most  easily 
calculated  by  setting  the  time  origin  ( t  —  0)  to  be  the  stopwatch’s  starting  time, 
then  determining  the  probability  that  there  are  no  impulses  in  the  time  span  [0,<]. 
It  follows  that  1  —  P(no  impulses  in  [0,  t])  is  the  probability  that  the  waiting  time 
is  less  than  t,  so  differentiating  this  quantity  results  in  the  pdf  fTw(t)- 

In  order  to  find  P(no  impulses  in  [0,f]),  the  pmf  of  No,  the  number  of  clusters 
active  at  time  zero,  must  be  calculated.  This  pmf  is  calculated  in  three  steps:  (i) 
The  number  of  clusters  N  created  in  the  range  [— T,  0]  is  a  Poisson(A2T)  random 
variable.  If  ordering  is  neglected,  the  starting  times  of  these  clusters  are  uniformly 
distributed  over  [— T,  0]  and  are  independent  of  each  other.  The  probability  Pa  that 
any  given  cluster  Xk  with  starting  time  4  ~  U([—  T,  0])  is  active  at  t  =  0  can  then 
be  calculated  by  conditioning  on  tk : 

Pa  =  ^(cluster  starting  in  [— T,  0]  is  active  at  t  =  0) 

=  f  J_T  PiTL  >  -r> dr 

-  <4'60> 

(ii)  Now  condition  on  N:  given  N  clusters  with  starting  times  uniformly  and  inde¬ 
pendently  distributed  over  [— T,  0],  each  with  a  probability  of  surviving  past  t  —  0 
given  by  (4.60),  then  the  number  of  clusters  surviving  at  t  =  0  is  a  binomial  random 
variable: 


P(k  clusters  active  at  t  =  0| AT) 

=  i  =  0,l,2...JV. 


(4.61) 
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(iii)  Since  N  is  Poisson^T),  the  final  result  is 


P(k  clusters  active  at  t  —  0) 


=  E 


(A  2r)”e 


n^—X  2T 


n\ 


Pj(l  -  PA)n~\  k  >  0 


n= k 

„-A2r 


n!  k\(n  —  k)\ 

PA  \*g(A  2T)k+n 

A '  n=0 


_  (  p*  V  v  IMXZn  _  p 

fc!  Vl-Pj  to  »!  1  ^ 


fc! 


(4.62) 

(4.63) 

(4.64) 


where  (4.63)  results  from  a  change  of  variables  n  =>  k+n.  From  (4.64)  it  is  clear  that 
No  is  a  Poisson  random  variable  a  parameter  equal  to  A2  T  scaled  by  the  probability 
(4.60),  a  well-known  result  for  a  Poisson  process  modulated  by  an  indicator  function. 
Setting  T->oo,  the  final  result  is  that  N0  is  a  Poisson(A2/A£,)  random  variable. 


Each  of  the  No  clusters  active  at  t  —  0  has  as  least  one  impulse  remaining.  In 
addition,  from  the  memory  less  property  of  exponential  random  variables,  the  first 
impulse  after  time  zero  in  each  of  these  clusters  has  a  time  of  occurrence  that  is 
EXP(Ai)  distributed.  Therefore,  in  order  for  there  to  be  no  impulses  in  [0,<],  all 
No  active  clusters  must  have  their  next  impulse  occur  after  t.  Conditioning  on  No 
results  in 


P(no  impulses  in  [0,  t\  due  to  clusters  starting  before  time  0) 


E 

71=0 

A- 


n\ 


(e-Y 


Taking  into  account  that  there  must  be  no  new  clusters  in  [0,  t]  as  well, 
P( no  impulses  in  [0,  £])  =  e~>'2te~JL  (1-e  1  ), 


(4.65) 

(4.66) 


(4.67) 


and  so  the  forward  recurrence  pdf  is  finally  found  by  negating  the  derivative  of 
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Eq.  (4.67): 


hwit)  = 


A2  + 


(4.68) 


4.4.2  Interarrival  Time  pdf 


The  interarrival  pdf  /y7(t)  calculation  is  similar  to  the  forward  recurrence  calculation 
except  that  there  is  by  definition  an  impulse  at  time  t  =  0.  It  is  not  known  whether 
this  impulse  is  the  first  impulse  of  a  new  cluster  or  an  additional  impulse  from  an 
established  cluster,  however,  so  it  is  necessary  to  condition  over  both  possibilities. 

First  calculate  the  expected  value  of  the  number  of  impulses  in  a  cluster,  which 
is  1  +  2£[iVfc].  Using  the  pmf  (4.14),  this  value  is 


1  +  £[AT*]  = 


oo 

i  +  E 


71—0 


nXL 

Ai  +  A  l 


1  + 


Ai. 

A  L 


Ai  +  Ax  / 


(4.69) 


Furthermore,  if  the  process  is  modeled  as  existing  for  all  time,  the  probability  that 
any  single  impulse  is  the  first  of  its  cluster  is  the  reciprocal  of  (4.69).  Using  both 
this  information  and  Eq.  (4.67),  in  addition  to  the  fact  that  a  new  cluster  at  time 
t  =  0  must  not  cause  an  impulse  in  (0,<],  results  in 


P( no  impulses  in  (0,  t]  given  an  impulse  at  t  =  0) 

- _ e~Nt  ,  -x2t  -42.(1— e-*i') 


\E[Nk  +  l]  E[Nk  + 1], 

Negating  the  derivative  of  (4.70)  then  gives  the  marginal  interarrival  pdf 
e-^-£(i-«-Alt) 


(4.70) 


hi{t)  = 


Ai  +  A  l 


•  ^AiA2  +  AjAx  +  A2A l  H — ^ —  c  ^lt  +  AiA2e  ,  t  >  0. 


(4.71) 
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Threshold 

Ai  (s'1) 

A2  (s'1) 

Ax  (s'1) 

5£[PC|] 

12.76 

4.71 

1.29 

10£[|X|) 

9.17 

2.84 

1.06 

15£[|X|] 

7.42 

2.09 

0.95 

20£[|X|] 

6.07 

1.71 

0.89 

Table  4.1:  Parameter  values  of  the  clustering  Poisson  pdf’s  in  Figure  4.1 


4.4.3  Data  Analysis 

Now  that  /tj {t)  has  been  determined,  it  is  compared  to  the  impulse  spacings  seen 
in  low-frequency  noise  data.  Figure  4.1  shows  the  fit  of  Eq.  (4.71)  to  the  interarrival 
time  pdf  of  Arrival  Heights  noise  data  in  the  25.5  -  27.5  kHz  range,  for  the  whole 
month  of  December  1995,  during  the  hours  08  -  16  UT  of  each  day.  The  four  lines 
in  the  plot  are  for  threshold  levels  of  5,  10,  15  and  20  times  the  expected  value  of  the 
magnitude  of  the  noise.  The  accurate  fit  of  Eq.  (4.71)  to  the  data  is  quite  apparent 
from  the  dotted  lines  in  the  figure;  the  parameter  values  used  for  each  dotted  line 
are  listed  in  Table  4.1. 

Many  more  data  samples  were  examined  in  addition  to  those  in  Figure  4.1,  but  all 
of  them  have  the  same  form  for  their  interarrival  pdf’s  independent  of  threshold  level 
(as  long  as  the  threshold  is  large  enough  to  distinguish  impulses  from  background 
noise).  The  term  A2  varies  with  location  and  seasonal  and  diurnal  variations,  but 
average  cluster  length  is  always  on  the  order  of  one  second  and  Ax  is  on  the  order 
of  5  -  10. 


4.4.4  Physical  Justification  of  the  Model 


Since  the  clustering  Poisson  equation  (4.59)  of  Section  4.4, 


y(t)  = 


00  Nk 

J2  2  *(*-** 


k=—oo  i— 0 


(4.72) 
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Figure  4.1:  Fit  of  the  clustering  Poisson  impulse  interarrival  pdf,  Eq.  (4.71)  (dotted 
lines),  to  Arrival  Heights  impulse  interarrival  data.  The  threshold  levels  are  5,  10, 
15  and  20  times  the  first  moment  of  the  noise. 


differs  from  equation  (4.4)  of  Section  4.2.2: 

Ne  NcNk 

',  =  E*  =  EE-7i%iM.  (4-73) 

k=l  k= 1  i=0  X 

it  remains  to  specify  the  connection  between  the  two.  Note  that  y(t )  in  (4.72)  is  a 
waveform  that  represents  the  detection  of  impulses  above  a  certain  threshold,  while 
Y  in  (4.73)  is  the  received  signal  at  a  given  point  in  time. 

Begin  by  defining  the  range  r  to  be  the  distance  from  a  given  cluster’s  source  to 
the  receiver  (r  =  |x|),  and  note  that  it  depends  only  on  |a|  and  r  whether  or  not 
an  impulse  crosses  a  given  threshold  Yth  in  the  received  waveform;  i.e.,  it  can  be 
specified  with  no  essential  loss  of  generality  that  cx  =  1  and  max  \E(rk,i]  &k,i)\  =  1, 
independent  of  6.  The  probability  that  an  impulse  from  the  cluster  crosses  Yth,  as 
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a  function  of  r  and  Yth,  is  then 

Prr*  =  P  (^T  >  «»)  =  2  Qr,  /«(“) da ■  (4-74) 

However,  this  poses  a  problem,  since  it  is  not  known  from  the  received  waveform 
how  many  impulses  and/or  entire  clusters  have  gone  undetected  because  of  small 
values  of  a.  In  order  to  proceed  further,  it  is  approximated  that  for  a  given  Yth, 
there  exists  some  cutoff  distance  rc  such  that  all  clusters  with  r  >  rc  can  be  ignored. 
This  cutoff  distance  is  chosen  to  be 


rc  = 


Ro 

Yh 


(4.75) 


for  some  constant  Rq. 

For  the  region  r  <  rc,  the  spatial  density  of  clusters  is  given  by  Eq.  (4.3): 

pM  =  (4-76) 

but  some  of  the  clusters  (especially  for  r  near  rc)  will  still  go  undetected  due  to 
propagation  attenuation.  To  compensate  for  propagation  effects,  an  effective  source 
density  is  defined  as 

pM  =  (4-77> 

Using  this  effective  source  density,  it  is  assumed  that  all  clusters  with  r  <  rc  are 
detected.  Thus  the  effective  number  of  spatial  cluster  sources  for  threshold  Yth,  in 
units  of  1/s,  is 

*- £**£-*'  <4-re> 

and  carrying  out  this  integration  gives 


A2  — 


PoCnRp  _  \ 

n  —  fj,  —  v  ) 


Ith  7 


(4.79) 


a  polynomial  in  Yth-  To  confirm  the  validity  of  this  approximation,  the  values  of 
relative  threshold  level  and  A2  listed  in  Table  4.1  were  fit  to  Eq.  (4.79)  by  taking 
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the  logarithm  of  both  sides  of  (4.79)  and  using  a  linear  regression.  The  exponent 
n  —  p  —  u  =  0.73  was  found  to  give  a  very  good  fit,  and  is  physically  justified  by  the 
expected  values  of  n  &  2,  u  «  1.0,  and  p  a  small  positive  number.  Thus  a  physical 
connection  between  the  parameters  of  (4.72)  and  (4.73)  is  shown.  (The  first  term 
in  parentheses  in  (4.79)  is  15.3,  but  this  number  is  irrelevant  since  po  is  unknown, 
Ro  is  arbitrary  and  only  relative  threshold  levels  are  used). 

It  is  seen  from  the  values  of  Ax  and  Ax,  in  Table  4.1  that  the  expected  number  of 
impulses  per  cluster,  1  +  Ax/A x,  is  on  the  order  of  5  -  10. 

4.5  Power  Spectral  Density  of  Clustered  Poisson 
Noise 

This  section  derives  the  power  spectral  density  (PSD)  of  clustered  Poisson  noise. 
It  is  shown  that  clustering  has  an  effect  on  the  power  spectral  density  only  when 
the  impulses  within  clusters  are  correlated,  so  for  the  clustering  Poisson  model  of 
Section  4.2.2,  the  PSD  is  simply  the  ensemble  average  of  the  PSD’s  of  individual 
impulses.  This  is  confirmed  in  the  data  by  calculating  PSD’s  of  (i)  isolated  impulses, 
and  (ii)  clustered  impulses,  and  showing  that  they  are  the  same. 

Consider  a  simple  Poisson  filtered-impulse  random  process 

OO 

r(t)=  £  Xk(t-tk),  (4.80) 

k=—oo 

where  Xk  is  the  random  waveform  of  the  kth.  event  and  tk  is  its  time  of  occurrence. 
Since  the  occurrence  of  events  is  a  Poisson  process,  the  random  variables  tk  —  tk- x 
are  i.i.d.  ~  EXP(A).  The  two-sided  power  spectral  density  of  Y(t)  is  given  in  [22] 
as 

SY(f)  =  A£[|*(/)|2]  +  (A£[*(0)]W),  /  €  K  (4.81) 

where  X(f)  is  the  Fourier  transform  (FT)  of  x(t), 

X(f)  =  r  x(t)e~i2*ft  dt , 

J  —  OO 


(4.82) 
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and  x(t)  is  a  given  realization  of  Xk(t)  (or  X(t)\  the  subscript  may  be  dropped  since 
the  Xk(tys  are  i.i.d.  waveforms). 


Now  define  the  X^s  of  Eq.  (4.80)  to  be  cluster  waveforms  similar  to  those  in 
Eq.  (4.4): 

Nk 

xt(*  -  «  =  E  E (‘  -  4  -  nr,  h,d,  (4.83) 

where  the  rfc>i  now  refer  to  the  position  of  the  ith  impulse  of  the  &th  cluster  relative 
to  the  start  of  that  cluster.  Then  let  Y ( t )  be  the  sum  of  all  the  clusters: 


k= 1 i=0  I**! 


(4.84) 


The  PSD  of  Y  is  still  Eq.  (4.81),  but  now  (4.83)  is  inserted  into  (4.80).  Note  once 
again  that  clusters  start  with  an  impulse  (the  i  =  0  terms),  and  the  Nk  s  are  i.i.d. 
geometric  random  variables  for  all  k. 


The  clusters  come  from  a  variety  of  locations  x  and  distances  |x|  and  they  are 
all  independent,  so  the  resulting  PSD  of  Y  is  the  summed  contributions  from  all 
locations.  Integrating  over  the  region  using  r  =  |x|  and  ignoring  the  zero  frequency 
(second)  term  of  Eq.  (4.81)  (it  is  assumed  to  be  zero)  gives 

Sy(f)  =  ^  dr^j  E  [|FT  of  any  one  cluster|2]  ,  (4.85) 

where  the  term  in  parentheses  is  the  effective  source  density  normalized  by  propaga¬ 
tion  attenuation,  as  in  Section  4.4.4.  (/?max  is  as  in  Section  4.3.2).  Further  defining 
this  term  in  parentheses  as  A  and  the  waveform  of  any  one  cluster  relative  to  its 
starting  point  as 

N 

X{t)  =  Yj  cia{E(t  -  n ;  6>{),  (4.86) 

2=0 

Eq.  (4.85)  can  be  written  as 


Sr(J)  =  A£[|  *(/)|2] 


(4.87) 
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roo  f 

\  .  .  I21 

A  E 

/  [Y^cia,iE(t -Ti]0i 

J~°°  \i= 0 

i 

CN 

1 

=  c\\E 
=  c\\E 


2=0 

£  £  «,•«,•£(/;  »<)£'(/; 

2=0 J=0 


(4.88) 

(4.89) 

(4.90) 


where  £(f ;  0)  is  the  Fourier  transform  of  F(f;  6)  and  (4.89)  results  from  the  Fourier 
shift  property.  In  order  to  proceed  further,  allow  the  a£s  to  be  wide  sense  stationary 
with  autocorrelation  function  R(k).  (They  are  still  mean  zero  with  variance  = 
72(0)).  Breaking  the  double  sum  into  terms  for  which  i  =  j  and  terms  for  which 
i  j  and  including  the  change  of  variable  k  =$■  |  j  —  i  |  in  the  latter  gives 


Sy(f)  =  Ac?(l  +  A)„lE  [| £(/;  9)|»]  +  cjA  I E  [£(/;  0)]|2 

~N—1  N-i 

*  £  J2  R(k)  (£[ei27r/(Ti+*-T,)]  +  F[e-j27r^(Ti+fc-T‘)]) 

.  2=0  k=l 


(4.91) 


However,  the  last  two  expectation  terms  are  the  characteristic  function  of  the  sum 
of  k  interarrival  times  within  a  cluster,  and  that  characteristic  function’s  complex 
conjugate,  using  2i rf  as  the  frequency  variable  instead  of  Since  the  interarrival 
times  are  i.i.d.  ~  EXP(Ai),  the  characteristic  function  expression  is  simply 


£[ej2ir/(T,+*-T0]  _ 


(4.92) 


where 

$(/)  =  r  Aie-AlV2^  dt  =  -  \  (4.93) 

Jo  Ai  -  j2irf 

Now  consider  two  different  cases  for  the  autocorrelation  Ra(k ):  For  the  first  case, 
set  Ra(k)  =  al5(k),  so  that 

W)  =  AcJ(l  +  A)^E  [|£(/;») |2 


(4.94) 
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In  this  case  the  PSD  of  the  signal  Y  is  simply  a  scaled  version  of  the  PSD  of  the 
individual  impulse  waveforms. 


For  the  second  case,  let  Ra(k )  =  cr2aftk\  where  0  <  /?  <  1.  The  double  sum  in 
(4.91)  can  be  converted  to  a  single  sum  such  that  Sy(f )  can  be  written 


SY(f)  =  Acf(l  +  ~)a^B  [|£(/;  0)|2]  +  cp  \E  [£(/;  0)]|2  a2a 


■En 


’  N 

J2(N  + 1  - *)  (0k$k(f)  +  Pk$k*{f)) 

.k= 1 


(4.95) 


which  indicates  that  $  is  simply  scaled  by  ft  for  this  case  of  R(k). 


In  order  to  solve  the  last  expectation  in  (4.95),  algebraic  summation  formulas 
are  used  to  find 


N 


£(1V  +  1  -  k)$k 

k—l 


(1V+1)(1  -$)$-$  + 
(1  -  $)2 


and  since  N  is  geometric  as  in  Eq.  (4.6),  it  then  follows  that 


En 


J2(N  +  1  - 


Uc=l 


(i  +  %)  (i  -  wmm-mf) + 
(i  -  mf)Y 


(4.96) 


(4.97) 


This  result,  when  inserted  into  (4.95),  is  cumbersome  to  solve,  but  it  can  be  simpli¬ 
fied  somewhat  with  the  approximation  /  >>  Ax,  which  is  true  for  the  VLF  frequency 
range.  (Note  from  Table  4.1  that  Ax  is  on  the  order  of  10).  Using  this  approximation, 


En 


E(iV  + 1  -  k)  (/?***(/)  +  /?*$**(/)) 

fc=l 

Ai  r  p Ax  _  p\i 


+ 


Al  LAi  -  j27r/  Ai  +  j2nf 


(4.98) 
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and  so 


Sr(f)  »  Acf(l  +  £ -)<rlE  [|£(/;0)|2]  +  |£[£(/; 


»)]f 


(4.99) 


Upon  examination  of  (4.99)  it  is  apparent  that  the  second  term,  the  one  due  to 
clustering,  is  negligible  relative  to  the  first  for  large  values  of  /.  Therefore,  in  the 
VLF  range,  clustering  has  no  effect  on  the  power  spectral  density  even  if  the  a,-’s  are 
highly  correlated.  To  verify  this  fact  in  the  data,  1/10  second  samples  containing 
one  isolated  sferic  (within  a  one  second  interval)  were  frequency  transformed,  then 
compared  with  the  frequency  transform  of  those  1/10  second  samples  containing 
many  sferics.  Averages  of  all  the  transforms  of  the  June  1986  data  at  Thule,  Green¬ 
land  are  shown  in  Figure  4.2;  it  is  clearly  seen  that  there  is  no  difference  in  the 
form  of  the  two  power  spectral  density  curves.  (The  y-axis  is  a  relative  index  only; 
the  plots  for  the  two  cases  are  offset  for  clarity.  The  rolloff  below  5  kHz  is  due  to 
a  high-pass  filter  inserted  to  reduce  power  line  harmonics).  Although  these  results 
do  not  provide  any  additional  proof  of  clustering,  they  do  show  that  clustering  does 
not  affect  the  PSD  of  a  signal  in  the  VLF  frequency  range. 


4.6  Conclusions 

This  chapter  proposes  a  statistical-physical  clustering  Poisson  model  for  atmospheric 
noise.  According  to  the  model,  the  sources  of  the  noise  are  clusters  of  impulses  with 
independent  and  identically  distributed  waveforms.  Cluster  occurrence  is  a  spatial 
and  temporal  Poisson  process  with  source  distribution  independent  of  direction  and 
time.  The  impulses  within  clusters  are  defined  as  variable  length  Poisson  processes. 
The  clustering  Poisson  model  is  justified  by  the  physical  properties  of  lightning,  and 
in  addition,  it  is  derived  in  conjunction  with  a  statistical  analysis  of  thousands  of 
hours  of  globally  collected  ELF /VLF /LF  radio  noise  data. 

The  model  is  verified  by  comparing  it  to  four  statistical  features  of  atmospheric 
radio  noise:  correlations  between  impulses  (Section  2.5),  envelope  amplitude  proba¬ 
bility  distributions  (Section  4.3),  impulse  spacings  (Section  4.4),  and  power  spectral 
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Figure  4.2:  Spectral  content  of  individual  sferics  versus  clusters  of  sferics.  The  data 
are  from  June  1986,  at  Thule,  Greenland.  There  is  a  lowpass  filter  inserted  near  7 
kHz.  The  spikes  below  5  kHz  are  power  line  harmonics,  and  the  spikes  above  10 
kHz  axe  other  man-made  signals  as  in  Figures  1.2  and  1.3. 
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densities  (Section  4.5).  The  correlation  properties  justify  the  assumption  that  the 
a's  of  Eq.  (4.4)  are  i.i.d.  and  that  the  clusters  and  impulses  within  them  are  in¬ 
dependent  as  well.  The  envelope  pdf  and  impulse  spacing  analyses  clearly  indicate 
that  the  theory  of  the  model  closely  matches  the  data,  and  finally,  clustering  of 
pulses  is  shown  to  have  little  effect  on  the  power  spectral  density  of  the  noise. 

Given  the  accuracy  to  which  the  predicted  statistical  features  fit  the  actual  data, 
the  clustering  Poisson  model  proves  to  be  a  good  candidate  as  a  model  for  atmo¬ 
spheric  noise. 


Chapter  5 

Low-Frequency  Communications 


5.1  Introduction 

As  discussed  in  Chapter  1,  there  are  a  number  of  digital  communication  systems 
that  operate  with  carrier  frequencies  in  the  VLF  and  LF  frequency  bands,  and  they 
employ  receivers  specifically  designed  for  impulsive  noise.  Since  the  best  receivers  are 
based  on  knowledge  of  the  statistics  of  the  additive  interfering  noise,  and  since  the 
previous  chapters  present  results  on  the  statistical  analysis  and  modeling  of  VLF/LF 
noise,  it  is  possible  that  these  results  may  be  used  to  design  better  performing 
receivers  than  those  currently  in  use. 

Specifically,  it  is  shown  in  Chapter  3  that  low-frequency  noise  has  an  a-stable 
envelope  distribution  at  locations  relatively  distant  from  heavy  sferic  activity.  Since 
many  VLF/LF  communication  systems  are  designed  to  operate  at  such  locations, 
it  is  of  practical  use  to  determine  whether  or  not  a  receiver  designed  specifically  for 
a-stable  noise  is  worth  its  complexity  in  terms  of  reduced  bit  error  rate  (BER). 

This  chapter  presents  receivers  designed  for  operation  in  a-stable  noise  and 
compares  their  BER  performance  to  the  conventional  nonlinear  analog  receivers 
commonly  used  with  impulsive  noise,  such  as  clippers,  blankers,  limiters  and  log- 
correlators.  The  signal  formats  examined  are  quadrature  phase-shift  keying  (QPSK) 
and  16  point  quadrature  amplitude  modulation  (16  QAM),  i.e.,  two  or  four  constel¬ 
lation  points  per  signal  dimension.  Hundreds  of  simulations  using  time-series  data 
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from  various  times  and  locations  and  at  various  center  frequencies  and  bandwidths 
are  performed,  and  the  following  results  are  found  uniformly:  for  QPSK  signals, 
virtually  no  performance  improvement  is  gained  when  using  an  a-stable  receiver 
instead  of  a  clipper  (the  best  conventional  receiver),  but  for  16  QAM  signals,  an 
improvement  of  several  dB  is  gained  by  using  an  a-stable  receiver.  In  addition,  it 
is  found  in  some  cases  that,  due  to  the  nature  of  impulsive  noise,  16  QAM  offers  a 
lower  BER  than  QPSK  for  the  same  average  transmitted  power  and  data  rate. 


5.2  Receiver  Structures  for  Impulsive  Noise 

This  section  assumes  basic  knowledge  of  digital  communications,  and  begins  with 
the  problem  already  cast  into  a  vector  form.  Reviews  of  digital  communication 
concepts  can  be  found  in  [32,  38,  43]. 

First  start  with  binary  phase-shift  keying  (BPSK).  The  receiver  must  choose 
between  two  equally  probable  hypotheses  h0  and  hi  (representing  either  a  zero  or  a 
one  bit)  based  on  the  received  vector  X: 

hi  :  X  =  S  +  N, 
ho:  X  =  — S  +  N, 

where  X,  N  and  S  are  complex  n-vectors.  This  problem  may  be  readily  generalized 
to  QPSK  in  two  dimensions,  noting  that  the  simulations  do  not  not  perform  symbol 
detection  ( i.e the  inphase  and  quadrature  channels  are  decoded  separately). 

For  4PAM  (pulse  amplitude  modulation),  the  problem  is  that  of  choosing  be¬ 
tween  the  four  choices 

hifi  :  X  =  3S  +  N, 
h1A:  X  =  S  +  N, 
h0A:  X  =  — S  +  N, 
hofi  :  X  =  — 3S  +  N, 

where  the  subscripts  on  h  indicate  the  bit  assignments  of  the  symbols.  The  gen¬ 
eralization  to  two  dimensions  results  in  a  16  QAM  constellation.  The  vector  X  is 
determined  by  sampling  the  analog  input  to  the  receiver  and  has  length  n  such  that 
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the  entire  nonzero  portion  of  the  original  symbol  s(t)  is  included.  If  this  nonzero 
portion  of  the  symbol  covers  time  T  (the  symbol  length)  and  the  receiver’s  bandpass 
front-end  input  is  limited  to  bandwidth  B,  n  is  on  the  order  of  n  =  2 BT.  Note  that 
B  refers  to  the  entire  front-end  detection  bandwidth  of  the  receiver,  which  can  be 
much  larger  than  the  bandwidth  occupied  by  s(t).  When  the  signal  is  downcon- 
verted,  the  inphase  and  quadrature  components  are  limited  in  frequency  to  and 
each  has  BT  degrees  of  freedom,  a  well-known  result  from  the  Nyquist  theorem. 

In  most  signal  detection  analyses,  N  is  assumed  to  be  an  uncorrelated  Gaussian 
random  vector  (GRV)  with  each  element  distributed  AT,-  ~  r?(0,  a2)  for  all  i  (where  r] 
refers  to  the  normal  distribution).  This  assumption  has  lead  to  a  number  of  elegant 
results  in  information  theory  pertaining  to  optimal  receivers  and  channel  capacity, 
but  it  is  not  applicable  for  atmospheric  noise. 


5.2.1  Maximum-Likelihood  Detection 


Using  BPSK  for  simplicity  and  given  that  hypotheses  hx  and  ho  are  equally  probable, 
the  decision  rule  that  minimizes  the  probability  of  error  for  a  received  vector  X  is 
to  compare  the  a  posteriori  probabilities  P(X|/ii)  and  P(X|h0)  and  choose  the 
hypothesis  corresponding  to  the  larger  one,  which  means  choose  hx  if 

/n(X  -  S) 

/n(X+S)  - 

and  ho  otherwise. 

For  N  an  i.i.d.  GRV  this  results  in  choosing  hx  if 


n 

2=1 


27 T£T2 


_ilal*ji>2 
e  2^ 


>1, 


27T<72 


which  can  easily  be  shown  by  taking  logarithms  and  rearranging  terms  to  be  equiv¬ 
alent  to  a  linear  matched  filter  detector,  i.e.  choose  h\  iff  (if  and  only  if)  X'S  >  0. 
The  matched  filter’s  impulse  response  is  the  time  inverted  signal,  demonstrating  the 
important  and  well-known  result  that  the  optimal  receiver  for  white  Gaussian  noise 
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(WGN)  filters  out  all  the  input  noise  except  for  the  noise  at  frequencies  occupied 
by  the  signal.  If  the  bandpass  signal  occupies  a  bandwidth  Bs  «  B,  then  the 
equivalent  baseband  signal  essentially  has  only  BST  free  dimensions,  and  only  the 
noise  samples  in  these  dimensions  affect  the  decision.  Two  important  results  that 
follow  for  WGN  are  that  (1)  the  front-end  bandwidth  of  the  receiver  need  only  be 
large  enough  to  pass  the  signal,  and  (2)  noise  and/or  other  man-made  signals  out¬ 
side  the  frequency  band  occupied  by  the  signal  of  interest  do  not  degrade  decision 
performance  because  they  are  rejected  by  the  matched  filter. 

For  impulsive  noise,  however,  the  noise  samples  are  not  Gaussian  and  noise 
outside  the  signal  bandwidth  must  be  considered  in  the  detection  process  if  possible. 
A  large  predetection  bandwidth  relative  to  the  signal  bandwidth  typically  enhances 
the  performance  of  VLF/LF  receivers  since  the  impulses  remain  sharp  relative  to 
the  symbol  length  and  can  be  more  cleanly  excised  or  compensated  for.  Otherwise 
the  impulses  are  spread  out  over  a  large  portion  of  the  symbol  waveform.  On  the 
other  hand,  the  predetection  bandwidth  must  be  small  enough  to  avoid  adjacent 
man-made  signals,  and  the  signal  bandwidth  must  be  large  enough  to  provide  a 
useful  information  rate. 

Once  the  signal  characteristics  and  noise  bandwidth  are  specified,  the  terms 
/n(X  —  S)  and  /n(X  +  S)  must  be  determined  in  order  to  implement  the  maximum- 
likelihood  (M-L)  receiver.  These  are  difficult  to  solve  for  impulsive  noise  distribu¬ 
tions,  however,  and  several  approximations  must  be  made.  Receiver  structures 
based  on  these  approximations  are  evaluated  in  Section  5.4,  but  first  a  review  of 
conventional  impulsive  noise  receivers  is  presented. 


5.2.2  Conventional  Nonlinear  Receivers 

Conventional  digital  receivers  for  VLF  noise  have  an  ad  hoc  zero  memory  nonlin¬ 
earity  followed  by  a  matched  filter.  The  nonlinearity  is  generally  either  a  clipper, 
blanker  or  hard  limiter,  although  simulations  show  that  the  clipping  amplifier  is  su¬ 
perior  to  the  other  two.  The  chosen  nonlinearity  operates  on  the  signal  plus  noise  in 
an  attempt  to  suppress  large,  short  duration  impulses,  thereby  making  an  attempt 
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to  “Gaussianize”  the  noise.  Of  course,  the  noise  must  be  of  significantly  greater 
bandwidth  than  the  signal  for  the  process  to  be  useful.  Since  a  matched  filter  fol¬ 
lows,  the  goal  is  to  eliminate  as  much  of  the  impulse  as  possible  from  the  band 
occupied  by  the  signal  while  minimizing  signal  distortion,  but  this  is  largely  a  trial 
and  error  process  since  both  the  receiver’s  front-end  bandwidth  and  the  nonlinear 
device’s  threshold  level  must  be  chosen  as  part  of  the  design.  Since  the  noise  is 
nonstationary,  one  or  both  of  these  parameters  may  have  to  be  adjusted  adaptively. 

The  chosen  nonlinearity  may  operate  on  a  wideband  receiver  input,  but  simula¬ 
tions  indicate  that  such  receivers  do  not  perform  well.  Although  many  sferics  stand 
above  all  the  man-made  signals  plus  background  noise,  the  effect  of  the  nonlinearity 
on  unwanted  man-made  interference  is  harmful  to  detection  (after  the  nonlinearity) 
in  the  band  of  the  signal.  For  this  reason,  the  predetection  bandwidth  is  chosen 
small  enough  to  eliminate  all  adjacent  man-made  interference,  and  the  nonlinearity 
is  applied  at  baseband  after  downconversion. 

Clippers  A  clipper  operates  linearly  unless  the  noise  exceeds  the  range  of  +/— 
a  given  threshold,  in  which  case  the  output  is  limited  to  the  +/—  threshold  level 
regardless  of  the  input.  Not  counting  signal  distortion,  using  a  clipper  results  in  a 
noise  pdf  that  is  truncated  to  the  threshold  value,  with  all  the  original  probability 
outside  of  the  threshold  range  now  placed  at  the  threshold  values.  The  matched  filter 
unit  that  follows  thus  sees  noise  that  is  somewhat  similar  to  a  Gaussian  distribution 
with  its  tails  cut  off. 

Blankers  Blankers  (or  hole  punchers)  simply  set  values  outside  some  threshold 
range  to  zero,  and  so  only  those  samples  uncorrupted  by  large  sferics  enter  the 
matched  filter.  Blankers  perform  poorly  when  compared  to  clippers  because  they 
discard  data  that  carry  significant  information,  even  if  those  data  are  corrupted  by 
impulses. 

Hard  Limiters  Limiters  output  one  value  for  inputs  greater  than  zero  and  the 
negative  of  that  value  for  inputs  less  than  zero.  They  are  simple  but  do  not  perform 
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as  well  as  clippers. 

Simulations  Of  the  three  conventional  nonlinear  receivers  just  described,  the  clip¬ 
ping  receiver  (a  clipper  followed  by  a  matched  filter)  consistently  performs  the  best. 
For  the  clipping  receiver  used  in  the  simulations,  a  number  of  clipping  threshold 
levels  are  used  but  only  the  one  giving  the  best  performance  is  shown  for  a  given 
plot. 

5.2.3  Previous  Receiver  Studies 

This  section  briefly  reviews  previous  work  in  developing  various  M-L  and  ad  hoc 
receivers,  all  of  which  perform  much  better  with  impulsive  noise  than  the  linear 
matched  filter  receiver  used  for  Gaussian  noise. 

Many  receivers  have  been  developed  using  a  small-signal  approach,  which  ex¬ 
pands  fx-s  in  a  Taylor  series  and  keeps  only  the  zero  and  first  order  terms  [30, 
39,  41].  Using  this  approximation,  the  M-L  receiver  is  comprised  of  a  zero-memory 
nonlinearity  followed  by  a  linear  detector,  and  so  it  is  similar  to  the  conventional  re¬ 
ceivers  described  above.  However,  the  nonlinearity  is  an  arbitrarily  complex  function 
that  generally  cannot  be  implemented  in  analog  hardware.  A  small  signal  approach 
is  often  inaccurate  because  many  low-frequency  communication  systems  have  quite 
large  signals  at  their  receivers  (see  Figure  4.2). 

Some  studies  choose  distributions  with  heavier  tails  than  those  of  the  Gaussian 
distribution  (e.g.,  Cauchy,  Weibull  or  beta),  determine  the  optimal  receiver  given 
the  distribution,  and  derive  receiver  results  either  empirically  or  with  computer  gen¬ 
erated  noise  [10,  31,  35].  Unfortunately,  the  distributions  used  are  often  inaccurate 
for  atmospheric  noise. 

The  one  M-L  receiver  with  possible  applications  to  this  work  is  the  log-correlator 
receiver,  the  optimal  receiver  for  noise  with  a  Hall  distribution  [23].  For  independent 
noise  samples,  the  optimal  receiver  structure  is  an  envelope  detector  followed  by  a 
logarithm  operation  followed  by  a  summation.  Since  the  Hall  envelope  pdf  is  found 
to  very  accurately  represent  impulsive  noise  under  certain  conditions,  it  is  expected 
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that  some  form  of  this  log-correlator  receiver  may  be  optimal  in  a  practical  VLF 
design.  For  this  reason,  the  log-correlator  is  included  in  the  simulation  analyses. 

Studies  of  actual  systems  using  a  clipper/matched  filter  combination  were  carried 
out  in  the  1970’s,  when  VLF  systems  were  first  being  installed.  Optimal  clipping 
levels  and  front-end  bandwidths  are  chosen  experimentally  and  are  system  and  lo¬ 
cation  dependent;  example  results  are  contained  in  [8,  11]. 


5.3  Alpha-Stable  Maximum  Likelihood  Receivers 

In  order  to  develop  an  optimal  M-L  receiver  for  a-stable  noise,  knowledge  of  the  joint 
distribution  of  the  elements  of  an  a-st able  random  vector  is  required.  Assuming  the 
elements  are  distributed  symmetrically  about  zero,  this  joint  pdf  /N(n)  is  defined 
by  the  joint  characteristic  function 

$N(t)  =  e/sl*'Bl0'‘<*),  (5.1) 

where  the  surface  S  is  the  unit  sphere  and  the  spectral  measure  fi(dS)  is  a  uniquely 
determined  finite  symmetric  measure  on  the  Borel  subsets  of  S  [2].  (The  variable  t 
is  equivalent  to  u_  in  Section  3.2.3.  Thus  while  a  mean  zero  Gaussian  random  vector 
(a  =  2)  of  length  N  is  uniquely  determined  by  an  N  x  N  covariance  matrix,  an 
a-stable  random  vector  ( a  ^  2)  for  even  N  =  2  requires  infinite  degrees  of  freedom 
to  be  specified.  Thus  it  is  impractical  to  determine  the  complete  description  of  an 
a-stable  random  process  from  the  noise. 

Even  though  several  simplified  types  of  a-stable  processes  exist  ( e.g the  sub- 
Gaussian,  linear,  and  harmonizable  forms  [2]),  none  are  applicable  to  the  noise  at 
hand.  In  order  to  proceed,  a  simplification  is  made  by  assuming  that  the  elements  of 
the  vector  are  independent  (but  not  identically  distributed),  an  assumption  justified 
by  the  statement  in  Section  2.6  that  no  appreciable  correlations  are  found  in  adjacent 
samples  of  the  background  noise.  Allowing  the  distribution  (i.e.,  the  parameters  a 
and  7)  to  vary  with  time  takes  into  account  the  rapid  statistical  variations  in  the 
noise  caused  by  bursts  of  impulses. 
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Covariation  As  an  additional  measure  of  the  lack  of  dependence  from  one  noise 
sample  to  the  next  in  VLF /LF  noise,  the  concept  of  covariation  is  introduced  and 
analyzed.  Second  order  statistics  such  as  covariance  functions  are  not  applicable 
for  a-stable  noise  since  the  variance  of  a-stable  noise  is  infinite,  so  the  covariation 
coefficient  is  used  instead.  It  is  specified  for  a  >  1  as  A xy,  where 

E(XY<>-'> ) 

A*y  -  mV)  (  * 

and  1  <  p  <  a  [40].  An  estimator  for  A  xy  using  p  —  1  is  then 


A  xy  ~ 


EL  Xisignffi) 

ZL  IK! 


(5.3) 


No  appreciable  covariations  are  found  for  time  differences  greater  than  1/B  (where 
B  is  the  bandwidth)  when  applying  (5.3)  to  the  data,  again  adding  validity  to  the 
assumption  that  adjacent  samples  can  be  considered  independent. 


5.3.1  Parameter  Estimations  of  Alpha-Stable  Distributions 

Since  the  noise  samples  are  assumed  to  be  independent  but  not  identically  dis¬ 
tributed,  it  remains  to  specify  the  two  parameters  a  and  7  of  the  a-stable  distribu¬ 
tion  as  a  function  of  time,  based  on  observation  of  the  noise.  Thus  it  is  necessary 
to  examine  parameter  estimation  techniques  for  a-stable  distributions. 

A  number  of  parameter  estimation  techniques  exist  for  a-stable  distributions  [37]; 
the  two  that  are  used  in  this  dissertation  are  (i)  histogram  fitting  and  (ii)  estimation 
of  the  characteristic  function.  The  former  is  explained  and  used  in  Chapter  3;  the 
latter  is  used  here  and  is  explained  as  follows:  the  characteristic  function  of  an 
a-stable  distribution  is 

(w)  =  E[ejwX]  =  e~7Ha ,  (5.4) 

and  thus  may  be  estimated  from  the  data  (assuming  a  symmetric  distribution  about 
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$x(w)  =  -77  ]Tcos(c<;Xi).  (5.5) 

iV  i= 1 

This  is  performed  in  the  analyses  for  ui  €  {0.5, 1.0, 1.5, 2.0}  and  then  a  linear  regres¬ 
sion  is  computed  using 


log(  log  $X  M)  =  log(7)  +  a  log  (a;)  (5.6) 

to  determine  a  and  7.  This  procedure  is  found  to  give  very  accurate  estimates  of 
a  and  7  when  compared  to  the  histogram  fitting  approach  (which  is  assumed  to  be 
the  most  accurate)  using  as  little  as  50  data  samples,  and  thus  is  used  throughout 
the  rest  of  this  work. 

Since  it  is  desired  to  estimate  a  and  7  adaptively  from  the  received  signal  plus 
noise  with  as  little  computation  as  possible,  two  other  parameter  estimation  schemes 
are  considered  as  well.  The  two  methods  examined  are  (iii)  to  derive  an  estimation 
based  on  the  number,  size  and  relative  times  of  localized  impulses,  and  (iv)  to  derive 
an  estimation  based  on  a  filtering  process  of  the  noise.  Many  possible  implemen¬ 
tations  of  method  (iii)  were  examined  but  none  accurately  predicted  the  localized 
values  of  a  and  7  determined  by  the  characteristic  function  estimation  method. 
With  method  (iv),  however,  a  simple  averaging  of  the  L  norm  is  found  to  correlate 
very  strongly  with  localized  values  of  a  and  7.  The  correlation  coefficient  of  the  lx 
norm  with  a  is  roughly  0.7  and  with  7  it  is  0.9. 

5.3.2  Realizations  of  Alpha-Stable  Receivers 

Four  realizations  of  a-stable  M-L  receivers  are  developed  and  simulated,  based  on 
the  parameters  estimation  techniques  of  the  previous  section.  There  are  endless 
receiver  specifications  possible,  but  the  four  used  here  are  typical  examples  and  are 
chosen  to  demonstrate  in  some  sense  the  best  and  worst  cases.  Of  primary  interest 
is  the  BER  improvement  realized  by  a-stable  receivers. 

Receiver  1  has  a  priori  knowledge  of  the  best  estimates  of  a  and  7  for  each  one 
minute  sample  of  data.  Changes  on  the  order  of  milliseconds  due  to  noise  bursts 
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are  not  accounted  for,  so  this  receiver  is  typically  the  poorest  performer  among  the 
four  a-stable  receivers. 

Receiver  2  has  a  p riori  knowledge  of  the  best  estimates  of  a  and  7  on  a  symbol-by- 
symbol  basis,  with  a  minimum  update  period  of  0.01  seconds  and  a  maximum  of  0.2 
seconds.  Short-term  estimates  of  the  noise  parameters  tend  to  improve  performance 
by  accounting  more  accurately  for  bursts  in  the  noise,  especially  for  faster  symbol 
rates. 

Receiver  3  estimates  the  parameters  a  and  7  on  a  simulation-sample  by  simulation- 
sample  basis  (the  simulation  sample  rate  is  2500/sec)  by  maintaining  a  100  sample 
running  average  of  the  /j  norm  of  the  noise  (without  signal).  This  running  average 
is  used  as  the  input  to  two  separate  linear  functions  to  determine  a  and  7;  the  gain 
and  offset  of  the  linear  equations  are  determined  a  priori. 

Receiver  4  is  identical  to  receiver  3  except  that  the  running  average  of  the  noise 
envelope  magnitude  is  calculated  from  the  signal  plus  noise.  The  noise  envelope  is 
estimated  by  subtracting  off  the  worst-case  (most  distant)  signal  point,  thus  giv¬ 
ing  an  indication  of  worst-case  performance  for  this  type  of  parameter  estimation 
scheme. 

The  four  receivers  described  above  are  somewhat  arbitrary,  but  they  are  demon¬ 
strative  of  the  performance  that  can  be  expected  from  practical  a-stable  receiver 
implementations.  Receiver  2  typically  performs  the  best  of  the  four,  since  it  has 
access  to  quick  and  accurate  estimates  of  the  noise.  Receivers  1  and  4  typically 
perform  the  worst,  1  being  slow  to  adjust  and  4  being  inaccurate  in  its  parameter 
estimation. 


5.4  Simulation  Results 

This  section  presents  the  results  of  hundreds  of  BER  simulations  performed  using 
actual  VLF/LF  noise  samples  as  additive  input  noise.  Signal  center  frequencies 
range  from  18  kHz  to  45  kHz  and  span  a  number  of  times,  locations  and  input 
bandwidths,  and  rectangular  pulses  are  used  as  the  symbol  basis  function.  The 
simulations  concentrate  on  those  times  and  locations  for  which  the  noise  is  most 
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accurately  modeled  as  a-stable,  but  other  cases  are  considered  as  well.  The  goals  of 
the  simulations  are  as  follows:  (i)  to  compare  clipping  receivers  to  limiters  and  linear 
(matched  filter)  receivers,  (ii)  to  show  the  effect  of  predetection  bandwidth  and/or 
symbol  rate  on  BER,  (iii)  to  compare  the  performance  of  a-stable  receivers  to  that  of 
clipping  receivers,  (iv)  to  determine  performance  differences  between  the  four  types 
of  a-stable  receivers  listed  in  Section  5.3.2,  and  (v)  to  compute  the  performance  loss 
when  using  16  QAM  rather  than  QPSK  with  a-stable  noise.  Note  that  (i)  and  (ii) 
confirm  previously  known  results. 

The  BER  curves  are  plotted  as  a  function  of  Eb/N0 ,  where  E\,  is  the  energy  per 
bit  and  No  is  twice  the  variance  of  the  measured  noise  samples  (so  that  No/2  is 
the  variance).  The  use  of  No  maintains  consistency  with  previous  literature,  but  is 
problematic  due  to  the  infinite  second  moment  of  a-stable  noise.  It  must  be  noted 
then  that  BER  comparisons  cannot  be  made  when  using  different  noise  samples; 
however,  BER  comparisons  using  the  same  noise  sample  (time,  station,  bandwidth 
and  center  frequency)  are  valid.  Both  two  and  four  level  constellations  are  examined 
(QPSK  and  16  QAM);  larger  constellations  are  not  considered  due  to  the  increased 
complexity  required  in  the  simulation  algorithms. 

Figure  5.1  shows  typical  simulation  results  for  QPSK  signals  in  impulsive  noise. 
The  noise  has  a  bandwidth  of  2  kHz,  centered  about  26.5  kHz,  and  is  from  Arrival 
Heights  during  June  of  1996.  The  symbol  rate  is  50/sec.  The  matched  filter  (the 
upper  solid  line)  is  clearly  the  worst  performer,  followed  by  the  limiter,  clipper,  log 
correlator  and  type  2  a-stable  receivers.  The  latter  three  are  very  close,  however,  and 
for  most  of  the  simulation  plots  they  are  essentially  identical.  This  result  confirms 
that  a  simple  clipping  amplifier  followed  by  a  matched  filter  is  nearly  optimal  for 
QPSK,  as  first  predicted  in  1974  by  Evans  and  Griffiths  [11]. 

Figure  5.2  shows  the  simulation  results  for  16  QAM  using  the  exact  same  noise 
samples  as  in  Figure  5.1.  In  this  case  the  two  M-L  based  receivers,  the  a-stable  and 
log-correlator  receivers,  are  clearly  optimal  to  the  conventional  nonlinear  receivers. 
Most  plots  show  the  a-stable  receiver  a  few  tenths  of  a  dB  ahead  of  the  log-correlator, 
and  the  latter  at  least  3  dB  better  than  the  rest  at  low  BER  (<  10-4). 

Figure  5.3  shows  the  performance  of  the  four  a-stable  receiver  types  described 
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Figure  5.1:  Typical  BER  curve  for  QPSK  noise,  comparing  linear,  limiter,  clip¬ 
per,  log  correlator  and  a-stable  receivers.  The  data  are  from  June  1996  at  Arrival 
Heights.  The  “linear”  solid  line  is  to  the  right  of  the  “stable”  solid  line. 


in  Section  5.3  for  a  16  QAM  constellation.  As  mentioned  above,  receiver  2  typically 
performs  the  best,  although  the  others  are  often  within  1.0  dB  at  BER’s  of  10-4. 
(The  2  dB  difference  seen  in  Figure  5.3  is  somewhat  of  a  worst  case). 

It  is  mentioned  earlier  that  decreasing  the  symbol  rate  or  increasing  the  predetec¬ 
tion  bandwidth  can  greatly  improve  performance,  especially  for  noise  characterized 
by  large  but  infrequent  impulses.  These  effects  are  shown  in  Figure  5.4,  where  the 
plots  show  BER  curves  for  symbol  times  from  1-100  symbols  per  second  and  the 
three  subfigures  themselves  show  variation  with  predetection  bandwidth.  (Note  that 
the  one  second  per  symbol  curve  on  Figure  5.4(a)  is  off  the  left  side  of  the  plot). 
Comparison  of  the  three  plots  demonstrates  the  large  variation  in  measured  No  that 
results  when  large  impulses  are  smoothed  by  narrower  and  narrower  predetection 
filters. 

The  final  topic  regarding  the  receiver  simulations  concerns  QPSK  versus  16 
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Figure  5.2:  Typical  BER  curve  for  16  QAM  noise,  comparing  linear,  limiter,  clip¬ 
per,  log  correlator  and  o-stable  receivers.  The  data  are  from  June  1996  at  Arrival 
Heights.  The  “linear”  solid  line  is  to  the  right  of  the  “stable”  solid  line. 


Figure  5.3:  BER  comparison  of  four  realizations  of  the  a-stable  receiver 
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QAM.  The  latter  requires  five  times  more  energy  than  the  former  in  order  to  main¬ 
tain  the  same  minimum  distance  in  the  signal  constellation,  which  translates  to  a 
5/2  (4  dB)  handicap  in  energy/bit  (as  clearly  seen  when  comparing  the  a-stable 
receivers  in  Figures  5.1  and  5.2  at  a  BER  of  10-5).  With  impulsive  noise,  however, 
performance  increases  with  symbol  time  for  the  same  Eb/N0,  so  it  is  possible  that 
this  gap  could  be  compensated  for  in  some  cases.  The  question  of  interest  then  is 
whether  it  more  advantageous  to  transmit  16  QAM  with  symbol  time  2T  or  QPSK 
with  symbol  time  T,  using  the  same  average  transmitted  power. 

Simulation  results  show  that  16  QAM  with  symbol  time  2T  does  in  fact  give 
better  performance  than  QPSK  with  symbol  time  T  for  high  Vd  noise  and  time- 
bandwidth  products  on  the  order  of  ten.  This  is  due  to  the  nature  of  impulsive 
noise,  which  when  compared  to  Gaussian  noise  with  the  same  No  will  have  a  much 
lower  level  of  background  noise.  Since  the  idea  of  signal  reception  in  sporadically 
impulsive  noise  is  to  evaluate  the  signal  “between  the  impulses”,  it  follows  that  a 
multilevel  constellation  could  be  transmitted  with  longer  symbols  to  take  advantage 
of  lower  level  background  noise,  with  no  loss  of  data  rate. 

Figure  5.5  shows  BER’s  for  three  cases  of  QPSK  versus  16  QAM.  For  low  and 
high  data  rates,  the  4  dB  performance  loss  is  clearly  seen,  but  at  125  symbols/second 
16  QAM,  the  curves  cross.  At  high  data  rates  (only  a  few  independent  noise  samples 
per  symbol),  a  large  impulse  dominates  the  symbol  regardless  of  the  signal  constel¬ 
lation.  For  low  data  rates,  the  greatest  percentage  of  data  samples  per  symbol  are 
background  noise  and  the  BER  curves  roll  off  in  a  shape  similar  to  those  for  Gaus¬ 
sian  noise.  For  BT  on  the  order  of  ten,  however,  the  low-level  background  noise  in 
a  double  length  symbol  can  make  up  for  a  short  impulse  somewhere  in  the  symbol, 
making  16  QAM  a  better  option  than  QPSK  for  the  same  transmitted  power. 


5.5  Conclusions 

This  chapter  shows  the  results  of  BER  receiver  simulations  performed  in  order  to 
compare  or-stable  receivers  to  common  analog  receivers.  The  conclusions  drawn  from 
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(c)  300  Hz  bandwidth 


Figure  5.4:  BER  curves  for  five  different  symbol  times  using  an  a-stable  receiver. 
The  noise  data  are  from  June  1986  at  Thule,  Greenland,  with  a  43.4  kHz  center 
frequency. 
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(a)  250  sym/s  16  QAM  vs.  500  sym/s  (b)  125  sym/s  16  QAM  vs.  250  sym/s 
QPSK  QPSK 


(c)  2.5  sym/s  16  QAM  vs.  5  sym/s  QPSK 


Figure  5.5:  QPSK  vs.  16  QAM,  with  16  QAM  at  half  the  symbol  rate  of  QPSK. 
The  noise  data  are  from  Thule,  Greenland,  with  a  43.4  kHz  center  frequency  and  a 
600  Hz  bandwidth. 
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these  simulations  are  as  follows:  (i)  For  VLF/LF  digital  communication  systems  us¬ 
ing  a  QPSK  signal  format,  an  a-stable  receiver  provides  essentially  no  performance 
improvement  over  a  conventional  clipping  receiver.  (This  result  can  likely  be  gen¬ 
eralized  for  all  other  constant  amplitude  constellations  as  well),  (ii)  For  a  16  QAM 
signal  format,  a  receiver  optimized  for  a-stable  noise  provides  several  dB  of  perfor¬ 
mance  improvement  at  low  BER.  (iii)  High  SNR  communication  systems  operating 
in  high  Vj  impulsive  noise  realize  1-2  dB  better  BER  performance  for  the  same  bit 
rate  and  average  transmitted  power  when  using  16  QAM  instead  of  QPSK. 


Chapter  6 
Conclusions 


This  chapter  summarizes  the  main  results  of  the  dissertation  and  discusses  possible 
areas  of  future  research. 


6.1  Summary  of  Results 

Chapter  2  introduces  some  of  the  statistical  properties  of  low-frequency  radio  noise. 
It  is  found  and/or  verified  that  (i)  seasonal  and  diurnal  variations  correlate  well  with 
global  lightning  patterns  and  can  be  used  to  track  global  climate  change,  (ii)  the 
narrowband  noise  envelope  pdf  of  low-frequency  noise  has  the  functional  form  of  a 
Rayleigh  pdf  at  low  values  of  dynamic  range  but  decays  with  an  inverse  power  law 
at  high  values,  (iii)  interarrival  time  distributions  indicate  clustering  of  sferics  with 
bursts  on  the  order  of  one  second  in  length,  (iv)  no  appreciable  correlations  are  found 
between  the  maximum  amplitudes  of  adjacent  sferics  or  between  these  maximum 
amplitudes  and  interarrival  times,  and  (v)  low-frequency  noise  can  be  viewed  as  the 
superposition  of  white  Gaussian  noise  and  discrete  noise  impulses  from  sferics. 

Chapter  3  describes  three  models  used  for  low-frequency  noise  envelope  pdf’s 
and  determines  the  parameters  and  accuracy  of  each  model  as  a  function  of  loca¬ 
tion,  time  and  frequency.  The  parameters  and  errors  of  each  model  are  found  to 
vary  with  thunderstorm  activity  and  the  noise  frequency  and  bandwidth.  The  chap¬ 
ter  concludes  that  the  Hall  model  is  the  optimal  choice  in  terms  of  accuracy  and 
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simplicity  for  locations  exposed  to  heavy  sferic  activity,  and  the  a-stable  model  is 
best  for  those  relatively  distant  from  heavy  sferic  activity. 

Chapter  4  proposes  a  statistical-physical  clustering  Poisson  model  for  atmo¬ 
spheric  noise.  According  to  the  model,  the  sources  of  the  noise  are  clusters  of  im¬ 
pulses  with  independent  and  identically  distributed  waveforms.  Cluster  occurrence 
is  a  spatial  and  temporal  Poisson  process  with  source  distribution  independent  of 
direction  and  time.  The  impulses  within  clusters  are  defined  as  variable  length  Pois¬ 
son  processes.  It  is  verified  that  the  clustering  Poisson  model  accurately  predicts 
the  statistical  features  of  the  actual  data,  and  thus  proves  to  be  a  good  candidate 
as  a  model  for  atmospheric  noise. 

Chapter  5  shows  the  results  of  BER  receiver  simulations  performed  in  order  to 
compare  a-stable  receivers  to  common  analog  receivers.  The  conclusions  drawn  from 
these  simulations  are  as  follows:  (i)  For  VLF/LF  digital  communication  systems  us¬ 
ing  a  QPSK  signal  format,  an  a-stable  receiver  provides  essentially  no  performance 
improvement  over  a  conventional  clipping  receiver,  (ii)  For  a  16  QAM  signal  for¬ 
mat,  a  receiver  optimized  for  a-stable  noise  provides  several  dB  of  performance 
improvement  at  low  BER.  (iii)  High  SNR  communication  systems  operating  in  high 
Vd  impulsive  noise  realize  1-2  dB  better  BER  performance  for  the  same  bit  rate  and 
average  transmitted  power  when  using  16  QAM  instead  of  QPSK. 


6.2  Topics  for  Future  Research 

Possible  topics  for  future  research  are  as  follows: 

•  Building  a  hardware  prototype  of  an  a-stable  receiver  and  using  it  in  an  ac¬ 
tual  communication  system  would  provide  useful  results  on  the  improvement 
realized  by  such  a  receiver.  The  receiver  would  likely  require  digital  signal 
processing  (DSP)  type  integrated  circuit  hardware. 

•  Obtaining  more  detailed  results  on  the  specification  of  the  a-stable  joint  pdf 
of  impulsive  noise  may  provide  opportunity  to  improve  a-stable  receiver  struc¬ 
tures  and  algorithms. 
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•  Acquiring  and  analyzing  a  continuous,  worldwide  collection  of  three  axis  ELF/ 
VLF /LF  time-series  data  with  extremely  accurate  system  timing  would  allow 
direction  finding  and  ranging  to  be  included  in  the  analysis  and  modeling.  Such 
a  database  will  be  collected  eventually,  and  will  allow  more  accurate  models 
of  low-frequency  atmospheric  noise  to  be  developed,  especially  regarding  the 
specification  of  cluster  waveforms  in  the  clustering  Poisson  model. 
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